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ABSTRACT 

When preplanetary bodies reach proportions of '^l km or larger in size, their ac- 
cretion rate is enhanced due to gravitational focusing (GF). We have developed a new 
numerical model to calculate the collisional evolution of the gravitationally-enhanced 
growth stage. The numerical model is novel as it attempts to preserve the individual 
particle nature of the bodies (like A^-body codes); yet it is statistical in nature since 
it must incorporate the very large number of planetesimals. We validate our approach 
against existing N-hody and statistical codes. Using the numerical model, we explore 
the characteristics of the runaway growth and the oligarchic growth accretion phases 
starting from an initial population of single planetesimal radius Rq. In models where 
the initial random velocity dispersion (as derived from their eccentricity) starts out 
below the escape speed of the planetesimal bodies, the system experiences runaway 
growth. We associate the initial runaway growth phase with increasing GF-factors for 
the largest body. We find that during the runaway growth phase the size distribution 
remains continuous but evolves into a power-law at the high mass end, consistent with 
previous studies. Furthermore, we find that the largest body accretes from all mass bins; 
a simple two component approximation is inapplicable during this stage. However, with 
growth the runaway body stirs up the random motions of the planetesimal population 
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from which it is accreting. Ultimately, this feedback stops the fast growth and the 
system passes into oligarchy, where competitor bodies from neighboring zones catch up 
in terms of mass. We identify the peak of GF with the transition between the runaway 
growth and oligarchy accretion stages. Compared to previous estimates, we find that 
the system leaves the runaway growth phase at a somewhat larger radius, especially at 
the outer disk. Furthermore, we assess the relevance of small, single-size fragments on 
the growth process. In classical models, where the initial velocity dispersion of bodies is 
small, these do not play a critical role during the runaway growth; however, in models 
that are characterized by large initial relative velocities due to external stirring of their 
random motions, a situation can emerge where fragments dominate the accretion, which 
could lead to a very fast growth. 

Subject headings: Asteroids — Origin, Solar system — Planetary formation — Plan- 
etesimals — Solar nebula 



Introduction 



One of today's key questions in planetary science is to understand the processes through which 
/im-size dust grains are converted into the ~10^ km-size bodies that constitute the rocky planets 
as well as the cores of the giant planets of our solar system. Nowadays this co llection has been 
enorm ously expanded with the discovery of hundreds of extrasolar planets (see lUdry and Santos 
20071 for a review). It is a challenge to understand a transformation process that spans more than 
35 orders of magnitude in mass, especially since observational data are rather limited. Although 
grain growth in protopla netary disks seern s to be a robust proce ss, as judged, for example, from 
(sub)millimeter studies (jNatta et al.l 120071 : iLommen et al.l l2009l ). the observational signature of 
macroscopic bodies becomes weaker and weaker with their size. In the solar system, important 
constraints include the remnants of objects that did not accrete into planets, like the bodies that 
constitute the Kuiper and Asteroid Belt as well as the meteoritic records that can be studied on 
Earth. 



In the core-accretion paradigm of pla.net f o rmation (|Mizuno et al.l 119781 : iHayashi et al.l Il985l : 



Lissauer and Steward 119931 : iPollack et al.l Il996l : iDominik et al.l 120071 ) two key ingredients play a 



crucial role in this transformation: molecular forces (surface forces) and gravity. The former 
dominates the behavior of dust particles at small scales, the latter at larger scales. At the high 
densities that characterize protoplanetary disks, micron size dust particles will qu ickly cluster into 
aggre gates - a process already observed in the cores of dense molecular clouds (jSteinacker et al. 
20ld ). Bodies of km-size have a gravity large enough to bind material and the same process 



allows these bodies to efficiently accrete each other. In this study we focus on the gravitationally- 
dominated regime. However, since it determines the initial conditions, we first briefly review 
theoretical efforts to overcome the problematic intermediate size range. 
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For meter-size boul ders surface forces seem to be too weak to act as an efficient sticking agent 
(|Blum and WurmI 120081 ) . Furthermore, the interaction with the (turbulent) gas especially affects 
meter-size particles. Whereas its damping properties are conducive to grow micron-size grains, the 
drag exerted on boulders causes them to drift inwards (and collide) at relati ve velocities of ^10 m 



s ^, except perhap s in some special settings like high-pressure environments (jKretke and Linll2007 



Brauer et al.ll2008l ). At these high velocities, laboratory dust collision experiments indicate these 
dust aggregates fragment rather than stick, or that perhaps growth already stops at the ~mm size 
(|Giittler et alJboiol : Iz^m et alJboiol ). 



Gravity has therefore been invoked to leap-frog the problematic meter-size regime. iGoldreich and Ward 



(jl973l ) argue that once the dispersion of dust particles drops below the threshold for gravitation- 
ally stability, the collapse produces ~km-size parent bodies. However, this mechanism ignores 
the role of gas drag, and of (self-induced) turbulence which will undoubtedly develop as a result 
of the angular velocity mismatch be tween the dust-dominated midplane layer and the pressure- 
supported gas (jWeidenschillingI Il980l ). In what is perhaps a surprising twist, recent studies have 
instead appealed to turbulence as a means to concentrate particles. For large, meter-size boul- 
ders the interplay with the gas leads to pile-ups of these particles (a.k.a. the streaming instability, 
Youdin and Johansenll2007l ) and numerical simulations have shown that these particles concentrate 
in cl umps, which subsequently collapse to form bod ies of ^10^ km to perhaps thousands of km in 



size (jjohansen et al.ll2007l . l2009l ). As another flavor. ICuzzi et al.l (j2008l ) argue that the intermittent 
properties of turbulence allow mm-size particles to concentrate in regions of overdensities of per- 
haps lO'^-lO^. Once captured by gravity, these clumps sl owly sediment to f orm the first generation 
of sandpile planetesimals, in the ~10-10^ km size range ( Cuzzi et al.ll2010l ). 



Thus, it is still unclear what the outcome of the primary accretion stage will be, i.e., the 
size distribution and the timescales. Recent work has tried to constrain planetesimal formation 
scenarios by comparing the outcome of models that cover t he gravitationally-domi nated phase 
with the preseii t day size distribution of the Asteroid Belt (jMorbidelli et al.l l2009l ). However, 
Morbidelli et al.l 1)20091 ) arrived at the conclusion that the Asteroid Belt population - in particular 
the 'bump' in the size distribution that is observed at at ~100 km (jBottke et al.ll2005l ) - is quite 
incompatible with accretion mod els, implying that it m ust have been the direct outcome of the 
primary accretion process (but see IWeidenschillingll2010l for a different interpretation). These works 
reveal the importance and motivation of modeling the gravity-dominated phase: its implications 
affect both the primary accretion process as well as the later stages of planet formation. 

After the formation of planetesimals growth becomes accelerated due to gravitational focusing 
(GF). GF increases the cross section for collisions by a factor ~(fesc/'u^)^ over the geometrical 
cross section, where w is the relative velocity between the bodies and Vesc the escape velocityEl For 
w < Vesc the accretion rate scales superlinearly with mass, dM/dt oc M ^, with k = 4/3, which causes 
bodies to separate in terms of mass over time (jKokubo and Idalll996l ). This phenomenon is better 



list of symbols and key abbreviations is provided in Appendix [Dl 
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known as runaway growth (jGreenberg et al.lll978l : IWetherill and Stewartlll989l ). The growth of the 
runaway body is relatively fast as long as the velocity dispersion v of the reservoir of (smaller) bodies 
from which it is accreting stays low, such that GF remains efficient, i.e., f ~ u; <C Vesc- However, 
GF does not only increase the collisional cross section, it also increases the rate of collisionless 
encounters, which dynamically heat the system, increasing w, and decreasing the focusing. 

Therefore, at some point runaway growth accretion will run out of steam. This transition has 
been stated by llda and Makind (| 19931 ) to take place at the point where 



2SmM > 



(1) 



with m and the mass and surface density of planetesimals and M and S m that of the runaway 
body(ies). Physically, this definition implies that the transition takes place at the point where 
the velocity dispersion of the planetesimal swarm is determined by the runaway body(ies). The 
consequence is that the growth rate (dM/dt) is reduced, although the growth mode stays locally 
runawayEI However, competitor bodies that are spatially separated are not affected by this heating; 
runaway growth in these zones is therefore faster than in the heated zones until these zones are 
also heated. The result is that (big) bodies that are dynamically separa ted tend to converge in 



terms of mass, a situation referred to as oligarchy (jKokubo and Idalll998l ). This dichotomy in the 
population - oligarchs (M) and 'field planetesimals' (m) - is a quite natural prediction of this 



for subsequent (semi-analytical) studies ( 


Thommes et al. 


2003; 


Goldreich et al. 


2006. 


2008; 


Fortier et al. 


2007; 


Brunini and Benvenutol 


200J 





To ensure sufficient growth within the lifetime of the gas disk (~10^ yr) the random motions 
of planetesimals v need to be damped during the oligarchy stage. When gas drag acts as the 
cooling mechanism it can be shown that the radius of the protoplanet Rpp increases only linearly 
with time, Rpp oc t, which for the ou ter disk means that the growth timescale becomes dangerousl y 



close to that of the gas disk lifetime (jKokubo and Idall2002l : iThommes et al 



2003 



Chambers 2006) 



Addi tional damping may be achieved through fragmentation ( ChamberslboOS : Kenvon and Bromley 



20091 ) or new, dy namically cold, reservoirs of planetesiinals c ould be tapped by the migration of 
the protoplanets (jTanaka and Idalll999l : iMordasini et al.ll2009l ). On the other hand, gap formation 
and/or gravit ational scattering could result in a strongly inhomogeneous disk in which accretion 
is suppressed (jRafikovl l2003l : iLevison et al.ll2010l ). Most of these models have used the assumption 
that the size distribution can be approximated by two populations (protoplanets and planetesimals). 
However, the validity of this assumptions (and others) depends ultimately on our understanding of 
the outcome of the runaway growth stage. 

Runaway growth (RG) calculations are quite challenging. The most straightforward way to 
assess the outcome of RG is by A^-body simulations. However, A^-body simulations, also suf- 
fer from severe computational constraints: A^ is restricted (typically to ~10^) and the dynamic 



^See Sect. 14.31 where we define runaway growth more precisely. 
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range (^M/m) these sim ulations can achieve is necessarily rather limited (|Kokubo and Idalll996l . 
2000; iBarnes et al.l 120091') . Clearly, to fo l low RG over a larger range one has to turn to statistical 
simulations (e.g., iGreenbe rg et alj Il978l . Il984l : IWeidenschilling et al.l 119971 : iKenvon and Luul 119981 : 



Inaba et al.l l200ll : iGlaschka i2006l ) . For a proper calculation at least four parameters should be 



followed: the mass m, semi-major axis a, eccentricity e, and inclination i. In most models, e and 
i are assumed to follow a distribution with parameters depending on mass and, if implemented, 
radial position only. The distribution function is altered by collisional (accretion, fragmentation) 
and coUisionless (viscous stirring, dynamical friction) processes and can be followed by numerical 
integration. A drawback of this approach is, however, that the distribution functions rely on large 
numbers: the number of particles of each type (m, a, e, i) should always be much greater than unity, 
whereas, as we saw above, in RG/oligarchy the distribution becomes discrete. In many codes the 
most massive bodies are therefore followed individually to incorporate th e 'local nature' of accretion 
(jWeidenschilling et al.lll997l : iBromlev and Kenvonll2006l : iGlaschkd l2006l ) . 



In this paper we describe a new method for the modeling of runaway and oligarchic growth. It 
is a Monte Carlo method in which each Monte Carlo 'particle' either describes a single body or an 
entire swarm of small bodies. We obviously start with the latter: when all bodies are still ~km size 
planetesimals each Monte Carlo particle represents a swarm of millions of planetesimals. However, 
the code preferentially favors bodies of larger mass; by the time a runaway body begins to dominate 
over the rest, it has already become a single Monte Carlo particle, and the individual particle nature 
of this runaway body is then automatically taken care of. In this way the transition from a fluid of 
planetesimals to a system with a few individual runaway bodies is properly treated. In fact, this 
would even allow a smooth transition to the next step of realism: an /^-body simulation, where 
the runaway bodies (oligarchs) would be treated as A^'-body particles. But we leave this next step 
to a later publication. 

With this new tool we embark on a parameter study of RG, in which the initial conditions, e.g., 
the mass and velocity distribution are varied. We include key physical processes like dynamical 
friction, gas drag, fragmentation, and turbulent stirring of bodies. Furthermore, we resolve the 
spatial dimension of the disk. Typically, simulations are followed until the biggest bodies reach ~10'^ 
km. We will define the timescale for runaway growth and identify the point where RG has been 
superseded by oligarchy. More generally, we assess the nature of gravitationally-dominated accretion 
under variation of the physical conditions. In particular, we address the role of fragmentation by 
adopting a very simple model where bodies colliding at velocities above their mutual escape speed 
convert a fraction of their mass into mm-size fragments. 

Section [2] presents the features of the multi-zone collision model. The collision model is vali- 
dated in Sect. [31 Readers more interested in the results may jump directly to Sect. [H which outlines 
the key characteristics of the runaway growth and oligarchy accretion phases. Section [5] presents our 
parameter study, where the accretion behavior is studied under variation of the physical conditions. 
Section [6] discusses some implications, while Sect. [7] summarizes the key results. 
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2. The coUisional evolution model 

In this section we address the key elements of the collision model. Section 12.11 provides some 
preliminaries. In Sect. [2^2] we describe our Monte-Carlo model, with which we calculate the time- 
evolution of the system. The zonal setup of the program is discussed in Sect. [231 Section [2.41 
introduces the interaction radii -Rjnt for the three velocity regimes under consideration. Section 12.51 
outlines how collisions are treated. Finally, Sect. 12.61 discusses merits and drawbacks of our ap- 
proach. 



2.1. Preliminaries 

For planetesimals - bodies which we consider to be of ~km-size or larger - we use the epicy- 
cle approximation to describe their orbital motion, in which the velocity is the vectorial sum of 
the Keplerian velocity Vk{a) corresponding to its semi-major axis a and a random component of 
magnitude v. When discussing interactions between two bodies of different mass we use vm for 
the random velocity of the largest body and Vm for that of the smallest. The orbital eccentricity 
is related to v as e ^ v/vk', similarly the random velocity in the vertical direction, Vz, is related to 
inclination i as Vz/vk — sini ~ i for i <C 1. 

The Hill sphere of a body of mass M and its corresponding Hill velocity are defined as 

/ M 

^h = af^] ; Vh = Rh^ (2) 

with the local orbital frequency and Mc the mass of the central star. The Hill radius approx- 
imately represents the distance over which 3-body interactions (the third body being the sun) 
become important. Using the escape velocity of the body, Uesc = 

^IGMjR, where G is Newton's 

constant, we find the useful auxiliary relation 

Rvl^^ = QRhvl, (3) 

which we will frequently employ. When discussing interactions between two bodies one should 
replace the above definitions by combined quantities, for which M = M1+M2 and R = Rg = R1+R2 
but this leaves Eq. [3] unaffected. 

Another usef ul dimensionless numb er is a, the ratio between the physical radius of a body and 



its Hill sphere fcf. iGoldreich et al.ll2004l ). 



Rh a \ M J ' \MqJ VAU/ 



where ps is the internal density of the body. Using Eq. [3] we find Vesc = f/i y^O/a. 
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2.2. The Monte Carlo collision model: interactions between representative bodies 

(RBs) 

2.2.1. Computational and physical particles 

Most statistical models involving planetesimal accretion rely on the concept of the distribution 
function to compute the collisional and dynamical evolution of a system. In this study, however, we 
calculate the evolution using a Monte Carlo method: a particle-based approach but still statistical 
in nature (rather than A^-body). At the core of the method is the concept of 'representative 
bodies' (RBs). These are the particles the computer program uses as a proxy for the full particle 
distribution of the physical system. One can say that each RB represents a group (or swarm) of 
physical particles. Because of their overwhelmingly large numbers, the behavior of each individual 
body cannot be followed; but with a limited representative sample - the RBs - a good census of 
the distribution can be obtained. 

The number of physical particles a RB represents is denoted Ng, the group size. There is 
complete freedom in choosing Ng] indeed, it will be a different number for each RB. In the computer 
program it is assumed that the physical particles corresponding to the RB share identical properties, 
i.e., the same masses, velocity dispersions, semi-major axis, etc. However, we do assume that the 
Ng particles are homogeneously distributed over the part of the simulation space the RB traverses; 
i.e., if the scaleheight of the RB is h^, the horizontal dispersion hx, and the semi-major axis a, we 
assume that the Ng physical bodies that are represented by the RB share the same o, h^, hz but that 
their phase angles characterizing the orbit are r a,ndomly distributed. It is impor t ant to point out 



that effects like reso nances or shepherding (e.g., IWeidenschilling and Davislll985l : [Patterson 



1987 



Mandell et al.l 120071 ) are not treated in our approach. 



Each representative body (RB) is characterized by four independent properties: (i) mass m; 
(ii) planar velocity dispersion (resulting from their eccentricity) v; (iii) vertical velocity dispersion 
(as resulting from their inclination) v^; (iv) semi-major axis a. Using these properties other particle 
properties are obtained, like the radius, R = (3m/47rps)^/^, or the escape velocity, Vesc- Each RB 
is further identified by a unique Ng. 



2. 2. 2. Choice and adaptation of Ng 

How is the relation between the physical particles and the RBs determined, i.e., the value of 
Ng7 'Traditional' MC-methods have Ng = 1: each computation body has a 1-1 correspondence 
to a physical body. This is, for example, sometimes used in aerosol coagulation studies (e.g., 



Gillespid Il975l ). However, for astrophysical purposes this is clearly inapplicable since then we 
only study the behavior of a limited number of physical bodies. If A'j.b.o is the initial number 
of RBs and rriQ the initial mass of the bodies, Ng = 1 implies that we only simulate a mass 
Mtot = NgN^y^fimo = N^y^fiTno, which, due to the modest amount of RBs computers can handle, 
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results only in a very limited potential for growth. Clearly, Ng ^> 1 is required, at least initially. 
In previous studies we h ave experimented with algorithms for choosing the proper group size 



Ng (jOrmel and Spaand l2008l ) . An algorithm in which the RBs represents a fraction of the total 



mass of the syste m (Mtnt), such that Nn oc Mtatim with m the rnass of the RB, is in most situations 
a good strategy (IZsom and Dullemondll2008l : lOrmel and Spaand l2008l ). However, it turns out that 
this mass sampling fails in systems that experience runaway growth. In such systems, a runaway 
body will form that separates from the continuous particle distribution. However, at th e point of 
separation the mass of the single run away body is typically negligible compared to Mtot (jWetherill 



199d : iMalvshkin and Goodmanll200ll ). For this reason, it is imperative that collision models resolve 
the high-mass tail of t he distribution very well. This 'unequal' mass sampling is the idea behind the 
'distribution method' (jOrmel and Spaanal2008l ). in which the RBs are distributed equally in terms 
of log mass, with the result that the high mass bodies are comparatively much better represented. 
In practice Ng = 1 is usually reached for the large bodies, whereas for the low- mass (~mo) particles 
Ng ^ 1. Perhaps counter intuitively, the total group mass of a low-m RB [Ngmo) is typically much 
larger than the mass of the high-m RB that has Ng = 1. 

In short, A'^ is a function of the RB's mass and since it is based on the current distribution, 
also a function of time. Furthermore, we recognize that in spatially isolated regions the distribution 
could evolve differently and let A'^ be a function of semi-major axis as well; hence, A^ = Ng{m, a, t). 
The way how the distribution method operates is discussed in Appendix O For a RB, we now 
distinguish between Ng, the number of bodies it represents at a certain point in the code, and A^, 
the number of representative bodies it should have according to the adopted algorithm for choosing 
Ng (here: the distribution method, Appendix lAj) . We desire that A^ = A^; however, N*{m,a,t) 
for a given RB varies over the course of the simulation run. How is the desired group size achieved? 

Let us illustrate. Suppose that the physical bodies associated with the RB undergo many 
accretion events such that they move up the mass hierarchy: m increases with respect to a charac- 
teristic mass of the system. The distribution method will then signify that they should be better 
resolved; A^* decreases and becomes less than A^, say A^* = Ng/2. The RB then splits into two 
identical RBs with Ng = N*, which from now on evolve differently. Likewise, if the RB is 'inert' 
(no accretion events) and declines in the mass-hierarchy, A^ will increase, perhaps becoming larger 
than Ng. We then say that the bodies are 'under resolved'. Such a situation is forbidden: allowing 
it would mean that the total number of representative bodies will increase indefinitely and strain 
computational resources. Another RB, sufficiently close in phase space, is sought to which the first 
RB is combined, averaging over their properties. However, this averaging is only done when the 
particle properties are indeed very close, v, Vz, and m should each be within 5% or less; otherwise 
the RB in question is kept in the program as a regular RB until such a situation does materialize. 

Due to these procedures, it is clear that the current number of RBs, A^rbi fluctuates with time. 
However, on average it is set by a target value - say, A"*^ - and this determines the resolution of 
the simulation. How A^ is determined form N*^ is explained in Appendix Rl 
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2.2.3. Interactions between the RBs 



How group encounters are implemented is described in detail in lOrmel and SpaansI (|2008l ) but 
let us illustrate the situation for a collision between two different RBs each representing Ngi and Ng2 
physical particles, respectively. The total physical collision rate is A12 = N giN g2/^Va,i2(^inc,i2 /y\at-, 
where /\va,i2 is the relative approach velocity, dine, 12 the interaction cross section, and l^nt the 
volume involved. Interactions among the Ng physical particles within a representative body are 
also accounted for and occur at a rate An = Ngi{Ngi — l)Afa^ii(7inc,ii/2Vint- The quantities 
At;a,c7int, and are all determined by the properties of the two RBs. Assuming, without loss 
of generality, that Ng2 > Ngi, the group collision is characterized by A'^int = Ng2 interactions in 
total with RB #1 accreting A^id^ = Ng2/Ngi physical particles of the second group. The collision 
then augments the mass of the first particle by an amount m2N\^^ and this becomes the new mass 
of RB #1 (mi — )• mi + A'idv"T'2)- The fact that typically A^int ^ 1 causes the code to speed up 
significantly. 

For m2 ^ mi and A^idv ^ 1 there is an important caveat, however. It is possible that Aidv'7i2 
is of the same order as mi, in which case we would instantaneously (i.e., during 1 group collision) 
increase the mass of RB ^1 by an amount of the order of its own mass, whereas in reality this 
occurs gradually. This is undesired since a sudden unphysical jump in the properties of RB #1 and 
in its collision rates with the other RBs is applied. Therefore, in the program we let the first (more 
massive) RB accrete a mass of at most /^mi and limit A^i^v accordingly, A^i^v ~ femi/m2. This 
means that only a fraction of the physical particles of RB #2 are involved in the group collision; 
the group is split. Note that the increase in the number of RBs by splitting up a single RB may 
cause Ng to fall below A^* (see above). The choice of reflects the computational cost. In this 
study, we adopt = 5 x 10^^ We have checked that our results are insensitive upon variation of 
fe by a factor of two. 



2.2.4. Summary of the collision code 

The flowchart of Fig. [1] summarizes the several steps of the program. The flowchart is intended 
to be schematic; it does not do justice to the full complexity of the underlying algorithm. First, the 
properties of each RB are defined, i.e., their masses, random velocities (inclination and eccentricity), 
and positions are assigned. We distinguish between collision rates and stirring rates, the latter 
determine the evolution of the random velocities. Steps B-D constitute the core of the program. 
Here, the collision partners are determined, together with the time step and the corresponding 
velocity changes. As explained below, in Sect. 12. 6^ we keep track of the random velocity change 
(Av^) of every RB by integrating its stirring rate {dv^ /dt) over time. After the collision has been 
performed in D, the particle properties have changed, which requires us to update the stirring and 
collision rates of all RBs in step E. The following steps, F-I, check several criteria. In F we check 
whether the RBs' random velocities have to be updated. If true, Av^ for the RB in question is put 
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(A) Initialization. For all RBs: 

- setup properties RBs 

(mass, velocities, position, etc) 

- calculate collision rates 

- calculate stirring rates dv2/dt 



z 



(B) MC step 




Using random numbers: 


- determine the 


two RBs 


that collide 




- determine timestep. At 






(C) Increment velocity changes 


For all RBs: 




Av2 = Av2 + 


(dv2/dt) At 



(D) Perform collision 

- determine number of 

colliders, Niciv 

- determine properties new RB 

- determine amount of fragments 

- determine new Ng of the RBs 




(E) update rates 



For all RBs: 

- collision rates 

- stirring rates 




Update Ng* — i 



Merge Rbs 

(if possible) 



Update velocity RB: 

v2 = v2 + Av2 
set Av2 = 



Fig. 1. — Flowchart describing the key steps of our colhsional evolution model. 



to zero. In G we check whether RBs have become under resolved, after which we merge RBs using 
the procedure described above. In H, finally, we check whether the function that determines the 
group size, N*, must be renewed. If any of these are positive, particle and stirring rates must again 
be calculated. After all theses criteria have been met, the state returns to B with a new cycle of 
the program. 



2.3. Spatial differentiation: the multi-zone setup 

Our code is multi-zone, with which we mean that RBs are assigned a particular semi-major 
axis on which their interaction radius depends. This situation is illustrated in Fig. [2j Here, we 
refer to the radial direction of the disk (the semi-major axis, a) as the x-direction and use z for 
the vertical direction. There are N^o zones each spanning a width Aa. The code then simulates a 
patch Osim = NzoAa in semi-major axis, centered at oq. The RBs are placed at the center of the 
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Aa 



Nnt 



1 



-4— <>- 



11 



h,=v,/Q 



,h^=v/Q 



Nzo 
— ► 



Fig. 2. — Sketch of the multi-zone setup. Symbols signify: Aa, resolution width; oq semi-major 
axis; Osim, simulation width; A^'^o, number of zones, hx,hz horizontal/ vertical scaleheights, -Rint, 
interaction radius. The scalewidth of the bodies is determined by the resolution of the grid Aa, 
hx = max{v/Q,, Aa). 



zones but (as mentioned before) assumed to be randomly distributed over the width of the zone, 
Aa, which sets the minimum value of the scalewidth, hx- However, if the horizontal excursions due 
to their eccentricities exceed Aa, like with the RBs in zone 11 of Fig. [2l their width is given by the 
eccentricities, or random velocities in the planar direction v, i.e., hx = v/Q. 

Whether or not RBs can interact (i.e., have non-zero rates) depends on whether their mutual 
excursions overlap (determined by hxi and hx2 and their radial distance) and on the interaction 
range, -Rint (discussed in Sect. 12. 4p . For example, if Rint = the RBs in zones 4 and 11 of Fig. [2]do 
not overlap and all rates are zero. On the other hand, if the RB in zone 4 is placed in zone 11 of 
Fig. [2] then they would perfectly overlap, and if it would be in zone 10 then there would be partial 
overlap. Also, if the interaction radius is large, say, i?int = lOAa, then there would be perfect 
overlap again, despite the fact that the RBs do not cross. In Appendix [C] the precise algorithm is 
presented. 

The way we treat the spatial dimension here echoes many f eatu res, albeit implemented more 
simply, of the multi-zone model developed by Spaute et al. ( 1991 ) and Weidenschilling et al. (jl997l ) 



(W97 in this section). But there are differences and we briefly mention these. In W97 when 
particles in bins are 'promoted' to individual bodies, these are assigned a true orbit (including 
phase angles), which is not done here. Another major difference is that W97 assume 'reflective' 
boundary condition, whereas in our treatment the boundaries are periodic: a body in the last zone 
Nzo hes computationally adjacent to the first zone. This means that in our code there is (for the 
moment) no radial gradient in the physical conditions, e.g., gas density, sound speed, etc. Our code 
is in that sense local. Like W97 when bodies merge they are placed in a new zone corresponding 
to the position of their common center of mass (since we treat a discrete grid, random numbers 
determine the zone the merged body is assigned to). However, for the remainder there is no spatial 
diffusion between the zones; e.g., there is no radial orbital decay of (small) particles due to gas drag 
or evolution in semi-major axis due to gravitational scattering. 
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Our implementation of the radial direction is (currently) not suited to model differences in 
evolution arising from global gradients in the physical conditions of the disk. However, even on 
scales where the physical conditions can be approximated to be constant, the disk will become 
inhomogeneous due to the emergence of runaway bodies, and it is this effect that we intend to 
explore in this study. Dynamical friction keeps these bodies rather cold and the low hx prevents 
them from overlapping. They become mutually isolated. A statistical 'particle-in-a-box' description 
does not do justice to this situation. Usually, single zone models resolve this problem by forcing the 



collision rates among the largest bodies to be zero (jWetherill and Stewartlll993l : llnaba et al.ll2001 



but this, perhaps somewhat ad-hoc prescription, is not adopted here. In addition, a runaway body, 
when sufficiently massive, starts to dynamically heat the planetesimal bodies, but only these in its 
neighborhood, which again makes the disk inhomogeneous (in phase space) lf| Thus, it is important 
to have a spatially resolved disk to assess the significance of these effects. 

How many zones are needed? Clearly, the more zones the better the spatial resolution and 
precision of the method. However, the number of zones will be limited by computational constraints 
because (i) each zone requires a minimum number of computational particles to resolve the mass 
spectrum and (ii) we compute the interaction between all RB pairs (not only these of the same 
zone). Using the properties of the system we define two key length scales. The first determines 
the resolution of the simulation (Aa) and is set to ^2t;esc,o/^) the spatial excursion the initial 
population of planetesimals (the bodies that contain most of the mass) would have if their random 
velocity equals the escape velocity, fesc,o- Here we anticipate that even if the initial random motions 
vq are <Cfcsc,o planetesimal-planetesimal stirring, which occurs on a timescale much shorter than 
accretion, will quickly heat up the planetesimals to velocities on the order of their escape velocities 
(jRafikovl I2OO3I ). However, if during runaway growth the dominant accre tion mode (in terms of 
mass) is between the biggest bodies (as postulated by lMakino et al.lll998l ) which are characterized 
by a very low hx we may still have oversampled the radial dimension. For these reasons we will 
test the dependence on resolution (Sect. 13. 3p . As the nominal resolution we use 



Aa 



2u 



osc,0 



SttGps R 



5 X 10"* AU 



Ps 



1 g cm 



-3 



1/2 



Rq 



10 km 



/ a \3/2 
VAU/ 



Note that this implies a much finer grid than used by lWeidenschilling et al.l (119971 ). where Aa 
AU is adopted. 



(5) 



0.01 



Similarly, we can set a condition for the total radial width of the simulation. Here, we anticipate 
that an olig arch or runaway body o f final mass Mj dominates a region several times its Hill sphere, 
e.g., 5-Rh,f (jKokubo and Idalll998l ). Such a length scale should fit comfortably within the total 
width of the simulation. Dividing the two lengths scales, 5-Rh,f and Aa, gives the minimum required 



^Here, 'neighborhood' is determined by the extent of the viscous stirring radius, R^b, and can exceed the width of 
the zone, Aa. 
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number of zones, A'^™, 

^min ^ 5i?h£ 1/6 / Ri/Ro \ -1/2 

where we used that Rh,i/Rh,o = Rf/Ro, M = AirpgR^/S and Eq. [5l In order to discriminate 
between an ohgarchic (several big bodies) or a runaway (one big body) outcome, N^o should be 
chosen several factors larger than A*"™™. From Eq. [6] it follows that for the same amount of growth, 
the inner disk requires more zones than the outer disk. 



2.4. The interactions: collisions, dynamical friction, and viscous stirring 

In our approach the collisional and dynamical evolution of the system follows from simple 
geometrical principles. For collisions this is a well-tested approach; the physical radii and relative 
velocity of two bodies directly determine its collision probability, or collision rate. We now extend 
this line of thinking and define an interaction cross-section (or radius, i?int) also for collisionless 
encounters. This may seem presumptuous since the gravitational interaction, being a long-range 
force, formally extends over an infinite distance; i.e., at any given time a planetesimal feels the 
force of many (strictly speaking: all) planetesimals. In this section we will only treat the close 
encounters that have the strongest influence on the orbit of a planetesimal, which result in a 
(finite) cross section for interaction. But for the flnal calculation of the stirring rates we will add 
a Coulomb factor to also include the more distant interactions (Appendix IB. 4p . We will introduce 
two cross sections (or interaction radii i?int) for the collisionless encounters: viscous stirring, i?vs) 
and dynamical friction, R^{. 

The interaction radii serve a twofold goal: (i) they determine the cross section of the interac- 
tions, which enter in the collision/encounter probability; and (ii) they determine whether bodies 
at different semi-major axes mutually influence each other. Every RB-pair is characterized by a 
unique R\^t and also a unique interaction outcome, e.g., the change in velocity is the same for all 
physical particles represented by the RB. 



The geometrical approach reflects iGoldreich et al.l (|2004l ) in their analytical study of oligarchic 



growth. It provides a very insightful treatment of how a population of planetesimals of mass m 
and a p opulation of olig a rchs of mass M with m <^ M mutually influence each other. However, 



whereas IGoldreich et al.l (|2004l ) treats a two components system, our model contains A^^b groups, 
amounting to ~-/Vr^b interactions, A^rb per RB. Our treatment is therefore numerical, but the under- 
lying principle is the same: the individual interactions between the groups follow from geometrical 
principles. 

We consider three different types of interactions: 



1. Collisions: accretion, bouncing, or fragmentation. The collisional radius is denoted, i2coi- 
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2. Dynamical friction. This occurs when particles experience a gravitational interaction but do 
not collide. This process leads to momentum exchange between the particles at an impact 
parameter b = R^f. In this work, i?df is given by the condition that the deflection angle 9 is 
large, ^ ~ 1. Then, the encounter can be approximated as an ID elastic collision. 

3. Viscous stirring. Apart from the exchange of momentum through dynamical friction, the 
nature of the encounter is such that energy can be extracted from or added to the Keplerian 
potential (see Appendix IB . 2 . 3P . This process is known as viscous stirring and operates at an 
interaction radius i?vs- The definition for the viscous stirring radius is set by the condition 
that the encounter changes the random velocity of the lightest particle Vm by a similar amount: 

AVm ~ Vm- 



The distinction between the dynamical friction and viscous stirring interactions should not 
be interpreted as meaning that these belong to two distinct encounters. In contrast, a (single) 
encounter will both contribute to the friction as well as the stirring. We simply dissect collisionless 
encounters into a part tha.t pres erves the random energy (dynamical friction) and a part that does 



not (viscous stirring) (jldalll990l ) 



For the calculation of the interaction radii we first need to specify the relative random velocity 
w (and Wz for the vertical direction), which is a function of the random velocities Vm and vm- 
Usually, the velocity of the lightest body Vm is larger than the velocity of the heaviest body vm 
due to dynamical friction but the situation may be different, e.g., in the presence of gas drag. 
The relative velocity at the point where the interaction takes place will in reality depend on the 
phase angles of the interacting bodies, and w will generally follow a distribution in velocity. These 
subtleties are ignored here and, as a crude approximation, w is simply taken equal to the maximum 
random velocity, i.e., 

w = max{vm;vM), (7) 
The relative random velocity w defines three velocity regimes: 

• the superescape regime, w > Vesc and Va = w; 

• the dispersion-dominated (d.d.) regime, 2.5vh < w < Vesc, and Va = w; 

• the shear-dominated (s.d.) regime, w < 2.5vh, and Va = 3b^}/2. 



Here, Va is the velocity at which the bodies approach each other. In the superescape and d.d.- 
regimes we have that the approach velocity, Va, equals w. However, if w becomes smaller than Vh, 
the s.d. -regime, the approach of the particles is determined by the Keplerian shear, Va ~ 350/2. 
Therefore, the gravitational regime {w < Vesc) splits into a d.d. -regime {vh ^ w < Vesc) and a 
s.d.-regime {w < Vh), see Fig. [3l Numerical studies have sh own that when w/z;;, <C 1 particles at 
impact parameters 6 = 2.5Rh will enter the Hil l sphere (e.g.. lNishidalll983l : iPetit and Henonlll986l : 



Ida and Nakazawalll989l : iGreenberg et al.lll99ll ). We adopt w = 2.5vh as the boundary separating 



the s.d.- and d.d. -regimes. 
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Fig. 3. — The interaction radii R^s, Rvs-d, Rdf, RcoU ^-iid the approach velocity Va as function of 
the relative velocity w between the interacting bodies. Radii are normalized to the Hill radius Rh 
and velocities to the Hill velocity Vh = Rh^ of the largest body. The dashed lines represent R^s 
and Va for distant interaction in the s.d. -regime. We have adopted a = Rs/Rh = 1-25 x 10"'^. 



Similar to Eq. [71 WG will dcfiiiG Wz — ixLajx('U^ 

,z',vm,z) as the relative velocity in the vertical 
direction. From this, we define the effective scaleheight h^s = Wz/^- Simplifying arguments like 
these are very common for the geometrical approach: they are perhaps not formally correct but 
for the moment they satisfy our purpose of a model that is accurate within factors of unity. Using 
the quantities for the interaction radii Rmt, the scaleheight hcs, and the velocity change upon 
interaction '^vf^^ we are able to construct interaction rates that only depend on these geometrical 
quantities. For example, for the stirring rates we obtain (see Eq. IB2p 

dv"^ TlRl^tRzVa 



dt 2/leff 



NsA^t, (8) 



where Rz is the interaction range in the vertical direction, Va the approach velocity (see below), 
and A'sj the column density of perturbers. Equation [8] can be compared with expressions that 
follow from more sophisticated studies. In Appendix [B] we quantify by how much our model is 
off, and adjust the collision rates that follow from our expressions accordingly, i.e., by inclusion of 
order-of-unity calibration constants into expressions like Eq. [HI 

We now provide expressions for Rcoi, Rdf, and R^g in these regimes. These are summarized in 
Table □ and Fig. [3l 
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2.4- 1- The superescape regime, w > Vesc 
Gravitational focusing (GF) is unimportant and all interactions are collisional iicoi = Rs 

Rl + i?2) Rdf — R-vs = 0- 



2.4-2. The dispersion- dominated regime, Vh ^ w < Vesc 

GF increases the collisional radius over the geometrical radius. The approach velocity is w 
but at impact the velocity is (at minimum) Vesc- Angular momentum conservation yields that the 
corresponding impact parameter at infinity is 

Rcol = Rs-^- (9) 
w 

Similarly, the criterion for dynamical friction is that the deflection angle changes over a large 
angle, or that /gAt ~ w, where fg ~ G{Mi + M2)/6^ is the gravitational force. Using At ~ b/w 
the corresponding impact parameter is therefore 

b=Rj!^y^R^,. (10) 



More formally, this impact parameter corresponds to a deflection angle oiw by ir/A (jBinnev and Tremaine 



20081). 



The reader could (correctly) argue that there is a certain level of arbitrariness in choosing i?df • 
For example, a deflection angle of = 7r/2 would amount to a dynamical friction radius that is only 
half that of Eq. \T0\ while for our formal deflnition of dynamical friction stated above - that it can be 
considered as a ID elastic collision - we would need 9 = n and b should be much smaller. However, 
Eq. [To] does give an indication of the scale at which strong interactions {9 ~ 1) become important. 
To complement the approach, as mentioned above, it is required to compute the resulting stirring 
rates {dv'^/dt, Eq. ED, compare these with existing literature treatments (Appendix [B]) . and, if 
necessary, to adjust the rate by the order-of-unity calibration constants. Somewhat surprisingly, 
the combination of Eq. [10] and the 'elastic ID collision' model (Appendix IB. 2p turns out to match 
very well the analytical result (Appendix IB. 4p . 

The criterion for viscous stirring is that the change in the random velocity of the lightest particle 
is significant, Avm ~ Vm- Thus, we solve fgb/w = Vm to obtain b = R^s = RsvlschmW, the radius 
for viscous stirring in the dispersion dominated regime. In the (usual) case that vm < Vm = w this 
equals the dynamical friction radius R^i, but if w = vm > Vm, Rvs > Rdi- 



2.4-3. The shear- dominated regime, w < v^; distant interactions 

In the s.d. -regime the approach velocity Va is given by the Keplerian shear instead of w. 
For these interactions the solar gravity cannot be neglected and the interaction includes three 
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bodies. Particles approachin g at distances b < will not enter the Hill sphere (jlda and Nakazawa 



IQsd : ICreenberg et al.l[l99ll l: instead, their trajectories strongly bend and the particles move away 



on horseshoe orbits. However, particles at slightly larger impact radii do enter the Hill sphere 
such that the characteristic impact radius is ^Rh- Following numerical and theoretical studies 



(e.g., iGreenberg et al.lll99ll ) we put the radius for entering the Hill sphere at 2.5Rfi- Similarly, 
the average approach velocity for particles entering the Hill sphere is calculated to be 3.2vh (see 
Appendix ETHl). 

Since particles at impact parameters b = 2.5Rh can be accreted, we put Rcoi = 2.5Rh- However, 
not every Hill-penetrating encounter results in a collision. In Appendix IB . 1 .31 we calculate the hit 
probability, /hit, with which expressions as the accretion rate must be supplemented. Alternatively, 
as a Oth-order approximation, we can define an 'effective impact parameter' by assuming the 2-body 
regime (Eq. [9|) and a relative velocity of ^2.5vh- Then, bco\ = R^oi^°'^^{w = 2.5vh) ~ \/RRh = 
a^^'^Rfi- This is what has been plotted in Fig. [3] but we emphasize that the program uses the -Rcoi~ 
/hit 'route', since this more accurately takes account of the spatial structure ( App endix IB . 1 . 3 p . 

Any particle that enters the Hill sphere experiences a strong interaction, such that i?df = 2.5Rh- 
For viscous stirring, on the other hand, significant stirring {Avm ~ fg'^'^ ~ ^m) already takes 
place at larger impact radii, R^s > Rh- Since the interaction timescale in this regime is set by 
the Keplerian shear, At ~ 17"^, the resulting impact parameter for these encounters becomes 
b ~ y/v^sc^/'^m^ = Rh\/Qvii/vm = Rm -d > Rh- These are lo ng-range forces that gravitationally 



perturb particles on non-crossing orbits (|Weidenschilline] Il989l ) 



Therefore, we distinguish between two viscous stirring radii in the s.d. -regime. When particles 
are capable to enter the Hill sphere, stirring is very efficient, because the velocity of these particles 
gets boosted to Vh, which can be 3>Um- For these particles R^s = 2.5Rh, Av"^ ~ vf^, and Va = 3.2vh- 
Otherwise, the viscous stirring radius is set to i?vs-d with an accompanying velocity change of (only) 
Av ~ Vm, and Va = 360/2 = SR^s-d^/'^ (see Fig. [3]). In the s.d.-regime we have that i?int ^» /leff 
and therefore Rz = h^s- Inserting these expressions into Eq. [8] we see that dv^/dt oc Rhv\ for 
Hill-penetrating encounters, while dv^/dt oc RhV^Vm for distant interactions. Thus, despite the 
fact that both Va and R^t are larger for the distant encounters, their heating rates {dv^ /dt) are 
less than the close. Hill penetrating encounters due to the 'boost' Avm ~ Vh the particles receive in 
the latter case. However, particles separated at impact parameters b > 2.5Rh can only be stirred 
by the distant interactions. 



2.5. The collision model, gas drag, and fragmentation 



In this study we adopt a collision model that contains key physical processes like accretion, 
fragmentation, and bouncing but is overall characterized by its simplicity. 
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2.5.1. Disk properties and gas drag 



Table [2] list the adopted disk properties. We will run simulations at three distinct disk radii. 
The mass of the central star is put at 1 Mq. The p arameters for the sound speed Cg and the nebula 
pressure parameter ry ~ (c n/vk)'^ are adopted from iNakagawa et al.l (Il986l l. following the minimum 
mass solar nebula profile (|Weidenschillinglll977bl : iHayashi et al.lll985l ). However, we v ary the solid 



surfac e density S and dust-to-gas ratio to enab le a comparison with the studies of llnaba et al 



(j200ll ) (for 1 AU) and iKenvon and Luul (jl998l ) (at 35 AU). Therefore, the underlying density 
structure does not strictly follow a power-law. Since we do not treat a global disk configuration, 
these deviations are not critical. 

In simulations with gas drag we apply a deceleration to the particle's velocity evolution on top 
of the accelerations that follow from gra vitational encount ers. We use the modified expressions of 
Adachi et al. ( 1976 ) as written down by ( Inaba et al.ll200ll ): 
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drag *drag 

where tdrag = SpsR/SCoPg^cL, is the particle friction time, E = 1.211, and Cd the drag coefficient. 
The drag coefficient equals Cd = 0.44 for large bodies but s mall bodies in the 35 A U simulations 
follow a different Stokes drag regime for which C^) is larger (IWeidenschillingi Il977al ) . 



2.5.2. Collisions and fragmentation behavior 



We assume that planetesimal bodies are rubble piles consisting of much smaller particles (frag- 
ments) of mm size. Due to the porous spaces, the internal density of the rubble piles is fixed at a 
low value of /o^ = 1 g cm^^, except of the mode ls at 1 AU where we have put ps = 3 g cm^^ to 



facilitate the comparison with llnaba et al.l (120011 ). Within a single (head-on) collision between two 
rubble piles the fragments will undergo many more collisions and dissipate a lot of the collision en- 
ergy, which renders the overall collision very inelastic. Thus, although th e individual coefficient o f 
restitution ei between two fragments is usually on the order of ~0.5 (e.g.. iHeiBelmann et al.ll2010l ). 
we assume that collisions between two rubble piles can be mo deled with a net e f fectiv e coefficient of 
restitution (e) that is much lower, e = 0.01 <^ 1. According to lGreenberg et al.l (119781 ) this value for 
the (effective) coefficient of restitution corresponds to loosely bounded regolith or weak material. 
In reality, e will further depend q i i the impact parameter and mor e sophisticated collision models 
are needed (jLeinhardt et al.ll2000l : iLeinhardt and Richardsonll2002l ). 



In our model, fragmentation only occurs when the relative velocity exceeds the escape velocity, 
Vesc- As a very simple prescription we assume the fragments are of the same size and very small. 



Table 1. Summary of the interaction radii i?int for collisional encounters (collisions) and 
collisionless encounters (dynamical friction and viscous stirring). 



Velocity regime Interaction radii, R[^t 

Collisions Dynamical friction Viscous stirring 

close distant 

i?col -Rdf Rvs Rvs-d 

Rs 

{6RsRh)'^\vh/w) 6Rh{vh/wf 6Rh{vl/v„,w) 

a^/^Rh 2.5i?h l.hRh i?h(6r;h/fm)i/2 



^Valid only when Rvs/w < 0, ^. Otherwise, the expression -Rvs-d for the shear-dominated 
regime applies 

'^This indicates the effective collision radius. The more general approach, which more accu- 
rately takes care of the spatial dimension, is presented in Appendix IB. 1.3[ 

Note. — w is the relative velocity in the d.d. -regime, Vm the random velocity of the smallest 
particle, and the Hill velocity of the largest particle. Equation [3] has been applied for the 
expressions in the d.d. -regime. 



Table 2. Adopted disk properties 



Parameter 


Symbol 


1 AU 


6 AU 


35 


AU 


Bodies' internal density 


Ps [ 


g cm~^] 


3.0 


1.0 


1.0 




Ratio R/Rh 


a 




5.2(-3) 


1.3(-3) 


2.1( 


-4) 


Solid surface density 


S [l 


5 cm^^] 


16.7 


2.0 


0.2 




Dust-to-gas ratio 






86 


56 


56 




Sound speed 


Cg [ 


cm s^"*"] 


1.0(5) 


6.2(4) 


4.1(4) 


Gas density 


Pz 1 


g cm~^] 


1.4(-9) 


9.5(-12) 


1.1( 


-13) 


Nebula pressure par. 


V 




1.8(-3) 


4.4(-3) 


1.1( 


-2) 


Drag coefficient 


Cd 




0.44 


0.44 


> 


.44 



W > Vesc 



Note. — Properties characterizing the physical conditions at 1, 6, and 35 AU. 
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Tf ^ Rq. Furthermore, we assume that the cohision dissipates the majority of the colhsion energy 
and that only an energy e-Ecol remains to eject the fragments, with Ecoi = mim2v'^/2{mi + 7712) 
the impact energjQ and e = 0.01. However, there is no energy required to break the material since 
the fragments are already loosely bound. Therefore, a mass fraction /^ag = eE'coi/l^T-tott'esc/^)' 
with mtot = ^1 + ^2 the combined mass of the collision partners, ends up as fragments with the 
mass of the main body being reduced correspondingly, M = mtot(l — /frag)- The choice of these 
parameters yields an impact strength for destruction of Qd — ^esc/2e ~ 10^ (i?/km)^ erg g~^, 
which f alls within the range of se veral proposed strength curves in this gravitationally-dominated 



regime (jBenz and Asphaug||l999l ). 



We will treat two values for the fragment size, afj. ~ 1 mm (chondrule size particles) and 
Ofrag = 10 cm (boulder- type particles). Fragments are not allowed to accrete among themselves 
but can be accreted by larger bodies. Next, we recognize that gas drag will influence the approach 
velocity of the fragments. We distinguish between two situations: (i) strong drag and (ii) weak 
drag. Strong drag occurs in the 1 and 6 AU simulations: the fragments are tied to the gas and move 
at a fixed relative velocity of v = r]Vk, corresponding to the subkeplerian gas velocity, and Vz = 0. 
This means in practice that the re-accretion of fragments is suppressed since their approach velocity 
Va is rather large with little or no GF. However, in the 35 AU models we relax the strong coupling 
assumption and model the dynamical behavior of the fragments in exactly the same way as the big 
bodies. Then, the gas drag only has a (slight) damping effect with the cooling being dominated 
by mutual (inelastic) collisions among the fragments. Since these collisions are abundant, these 
are usually very effective to dissipate any random motion, despite the stirring of the big bodies to 
which the fragments are also subject to. 



2.5.3. Turbulent stirring 



by turbulence 


(Lauehlin et al. 


2004 




Nelson 


2005 


). This results in an eccentricity evolution 


( Oeihara et al. 


2007; 


Ida et al. 


2008) 





e ~ O.I7 



^91 



a \ ' 



{'AUJ \Tk 



t 



1/2 



(12) 



where Sgi = 2400 g cm ^. Inserting this value for Sgi, Tk = 2Ti/VL for the orbital period, and 
squaring gives 

17, (13) 



dt 



1.6 X 10^^72 



90 Mff 



^Since we treat the Va > Vcsc case, we neglect the focusing term for the impact energy. 
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which in terms of velocity units v = eaVL reads 



dt 



3.5 X 10 



cm^ s ^ 



7 



10' 



100 Ma 



/ a \-5/2 
VAuJ 



(14) 



In the above expressions the dimensionless 7 embodi e s the strength of the turbulent density fluctu- 



ations. Using 7 in the ran ge of 10 ^-10 ^ Ida et al.l ( 20081 ) show that most planetesimal colhsions 
result in destruction (e.g., iMorbidelli et al.ll2009l ). These large values for 7 were based on global 



simulations invol ving the magneto- rotational instability (MRI). However, subsequent local shearing 
box simulations (jYang et al.l l2009l ) indicated a lower value, 7 ~ 10~^, despite the fact that these 
MRI simulations gave a rather high turbulent-alpha parameter of ot ~ 10-2. Therefore, we will 
adopt both 7 = 10"^ and 7 = 10""^ when including turbulent stirring in our simulations. However, 
provided an sufficient shielding by (sub)//m size grains, th e turbulence may be si gnificantly sup- 
pressed in the interior regions of the disks (dead zones) (e.g. 
stirring may perhaps not be so effective. 



Turner and Sandl2008l ) , and turbulent 



2.6. Merits and drawbacks of the collision model 

We end this section with an assessment of the collision model, discussing its strengths and 
weaknesses, and sketch avenues for future extensions. As the strengths of the collision model we 
list: 

1. the particle nature of the geometric model; 

2. the large dynamic range concerning runaway growth; 

3. the presence of stochastic effects due to the Monte Carlo noise; 

4. the intuitive nature of the geometric model. 

Here, the first three points are related to the numerical model (Sect. [2?2]) . whereas the last con- 
cerns the interactions (Sect. [2^ . We find the way we treat interactions, which is inspired by 



Goldreich et al.l (|2004l ) 2-group's approximation, more intuitive tha n the rather comp lex (but per 



3r compi 

haps more precise) formalisms of existing statistical programs (e.g. Jlnaba et aPboOll ). The advan- 
tage of the geometrical approach is that it identifies the critical mechanisms that drive the evolution 
in a transparent way - just because the key ingredients are all physically-intuitive properties like 
length scales and velocity changes. However, this is primarily a matter of taste; our numerical 
model would work just as well with the more formal expressions for the interaction rates. 

Arguably the biggest advantage of our Monte Carlo approach is that it deals with particles - the 
representative bodies (RBs) ~ which properties are independent of each other. For example, we can 
have two RBs with the same mass but with different velocities due to a sudden stochastic encounter. 
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This is quite different from the usual mass-binning methods, where the mass of the bin determines 
ah other properties. For the MC-method it is also easier to include more particle properties, for 
example, the internal structure of the bodies (molten or primordial) without any direct increase 
in the computational costs. Its particle nature and the many (independent) properties with which 
RBs can be quantified render the collision model especially attractive. 



There are many codes that mix statis tical and discrete elements (jWeidenschilling et al.l 119971 : 



Bromley and Kenvonll2006l : lGlaschkdl2006l ). Usually this involves a transition mass at which bodies 
are 'promoted' to an iV-body particle or a discrete particle. Our model, in terms of the representa- 
tive bodies, offers perhaps a more natural way to implement the transition; i.e., individual bodies 
are those which have Ng = 1 and this does depend on mass. Indeed, Ng = 1 bodies are formed 
pretty early in the course of the simulation run, a prerequisite to properly follow the runaway 
growth process. 

As the key drawbacks/omissions we list: 

1. the inefficiency of the Monte Carlo algorithm; 

2. the difficulty to model (semi) steady-state systems; 

3. the limited spatial diffusion of particles. 

The most severe disadvantage of MC-methods is that they are computationally very inefficient: 
of the ~iVj^b collision rates, only ~A'i.b are used. A run of our model, although much faster than 
A^-body, takes much longer than statistical m ethods based on the Smoluchowski (mass-binning) 



approach and is additionally rather noisy (see lOkuzumi et al.l |2009| for further discussion) . As a 



consequence, we were unable to explicitly calculate the velocity distribution within a swarm of 
bodies (i.e., to relax the assumption of a fixed distribution; see Eq. IB3P since this would have 
required too many RBs. 

For the same reason, it is too demanding to treat collisionless encounters also on an event- 
based approach. Due to their increased gravitational focusing (Fig. [3D collisionless encounters are 
more frequent than collisional interactions by several orders of magnitude, especially when the 
system relaxes to a quasi-steady state in velocity space. The problem is that the MC code does not 
recognize that the collective effect of these encounters cancels out, but instead resolves the strongly 
fluctuating velocities of the bodies, which render the code very inefficient. Therefore, the random 
velocities - eccentricities and inclinations - of the swarms are updated in a continuous fashion, as 
described in Fig. [TJ For every RB the cumulative effect of all Aj-b interactions is calculated, resulting 
in a stirring rate, dv'^/dt, and a (cumulative) velocity change Av'^{t). And only if the relative 
incremental change has exceeded its current value by a few percent (i.e., \Av'^\/v'^ = ~ 0.05) is 
the particle's velocity updated {v'^ v'^ + Av^), together with all ~Arb interactions it is involved 
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Thus, the velocity evolution of a swarm of bodies v{t) evolves smoothly with time. Although 
this is fine for the big bodies, strong scatterings are inherently discrete for small bodies. This 
especially concerns the s.d. -regime, in which the velocity of small bodies jumps to (which is 
possibly ^v) after a single scattering; i.e., its velocity-evolution is spiky rather than continuous. 
With our MC-model we have the tools to address these drawbacks, however, and in future work 
we may switch to an event-based evolution of the velocity (provided the associated increase in 
computational effort can be overcome). 

Although in the collision model we have resolved the radial direction of the bo dies, any spatial 
diffusion of planetesimals or ruii away bodies due to gravitational encounters (e.g., iTanaka and Ida 
19971 : lOhtsuki and Tanakal l2003l ) were ignored. Likewise, scattering by turbulent density fluctua- 
tions will further contribute to the radial mixing. These effects may render the distribution more 
homogeneous over large distances than our zonal setup currently supposes. On the other hand, if 
scattering by the oligarchs /runa way body dorninates, the surface density distribution of the disk 
may also become inhomogeneous (jRafikovll2003l : lLevison et al.ll20ld ). Gap formation is not included 
in the current setup of the model; but it could become important for the oligarchy stage. To study 
these and later stages the width of the simulation Aosim has to be increased, which invalidates the 
local assumption we apply (Aosim/o <^ 1). On a global (AU) scale, the differences in the evolution 
timescale may stifle the onset of run away growth in the outer disk, due to the effect of long-term 
perturbations (IWeidenschillingI l2008l ) . These effects cannot be self-consistently explored with our 
current setup. However, we will assess the consequence of a pre-stirred population of planetesimals 
by considering the superescape regime in which the (initial) random velocity of bodies exceed their 
escape velocity v > Vesc- 



Model validation 



3.1. Comparison with A^-body: viscous stirring and dynamical friction 



In this test, we copy the setup of the A^-body simulations of IStewart and Idal (|2000|) (SIOO 
in this section). In the SIOO A^-body simulations the dynamical behavior among planetesimals is 
studied without accretion. Therefore, we switch off accretion in our code (i?coi = 0) and only treat 
dynamical friction and viscous stirring. 



^We confirmed that the results were insensitive to a change of a factor 2 of the control parameter f^. 
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Fig. 4. — Test of viscous stirring. The evolution of the random velocity component (y-axis) of a 
population of 800 equal-size (m = lO^'* g) planetesimals as function of time (x-axis) for both a linear 
{left panel) and logarithmic {right panel) scaling. The planetesimals are distributed over an annulus 
of 0.056 AU centered at 1 AU. Eccentricity evolution is given by the black curve, inclinations by 
the grey curve. The dashed curve in panel B gives the ratio of inclinations to eccentricities, or 
/3 = Vz/v. The dotted horizontal line signifies the tr ansition between the s hear-dominated and 
dispersion-dominated regimes. Compare with Fig. 4 of Stewart and Ida (j200ol ). 



3.1.1. Equal mass system 

In the first test, 800 m = 10^*^ g planetesimals are positioned in a narrow annulus of 0.056 AU 
at a distance of 1 AU from the sun. Since bodies are of equal mass and no accretion takes place, 
the evolution of the system is determined by viscous stirring. The bodies initially have a very 
low random velocity, 1 cm s~^, but this increases with time due to viscous stirring. The resulting 
velocities are plotted in Fig. HI on a linear scale (panel A) and a logarithmic scale (panel B) and are 
normalized to the (combined) Hill velocity {v^ ~ 21 m s^^). Comparing Fig. with Fig. 4 of SIOO 
we find that the curves match to the 10% level (with most of the discrepancy in the inclination). 
This confirms the validity of our viscous stirring expressions, at least for the dispersion-dominated 
regime which we are probing with this test. 

Figure Hb shows the evolution resulting from our model over a more extended domain in time, 
which is instructive since it displays the underlying expressions of our geometrical model. In the 
s.d. -regime {v <C v^), it shows that eccentricities grow much faster than inclinations, causing the 
/3 parameter {/3 = Vz/v) to decrease. The reason for this behavior is that the interaction radius 
of viscous stirring i?vs is larger than the scaleheight of the planetesimals, R^s ~ -R/i ^ ^eff and 
that therefore the interactions take place in a planar (2D) geometry, which suppresses the vertical 
stirring. After i ~ 10 yr, however, the interactions enter the d.d. -regime, v > 2.5vfi. At this 
stage the interaction geometry becomes three dimensional, which means that f5 evolves towards its 
equilibrium value, /3 ~ 0.5. This effect initially stagnates the eccentricity stirring at the expense 



- 25 - 



0.015 



0.010 



0.005 



0.000 - 



m = 10 g bodies 




X 



X 







5000 
Time [yr] 



10000 



Fig. 5. — Like Fig. |3]but for a two component system of m = 10^^ g (grey curves) and m = 10^^ g 
bodies (black curves) with an equal amount of mass in both components. The tot al surface density 
i s aga in S = 10 g cm~^. Inclinations are given by dashed curves. Compare with IStewart and Ida 
(l200d ). their Fig. 9b. 



of the inclinations but after t ~ 10^ yr the equilibrium has been achieved and v and Vz evolve on 
similar timescales. Finally, after t ~ 10^ yr the interactions have reached the superescape regime 
{v = Vqsc) and viscous stirring does not operate anymore. Bodies cannot be stirred above their 
mutual escape velocitylf] 

Equation [8] can be used to understand the qualitative behavior of the curves. In the s.d.-regime, 
where the disk is thin, we have that Rz = /igfr, Rx ~ Rh, Va ~ Vh, and Au^ = v^. Since these 
quantities are all constant we find that the random velocities grow with the square root of time, 
V oc t^l"^. In the d.d. -regime, on the other hand we have that RxRz ~ R^{vesc/vY ■, Va/heS constant, 
and Av"^ = v"^ . This results in dv'^/dt oc v~'^ and therefore v oc These slopes are indeed 



obser ved in Fig. UJa and agree with more detailed previous studies (e.g., Ilda 



199C 



Ida and Making 



3.1.2. Two component system 

Next, we calculate the evolution of a two component system, adopting the same parameters as 
before, but replacing half of the mass of the 10^^ g bodies with bodies of 10^^ g. Thus, we have a 
two-component system in which we expect dynamical friction to operate. Our results are displayed 
in Fig. [5] and should be compared with Fig. 9b of SIOO. Again, we find a good match at the 10% 



®In this test we have fixed the internal density of the bodies at = 3 g cm However, if these are truly point 
sources (infinite ps) the superescape regime does not exist and the flattening would take place. 
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Fig. 6. — The cumulative number distribution and velocity distr i butio n at several times during 
the evolution at 1 AU. Parameters are the same as in llnaba et al.l (120011). compare with th eir Fig. 
9. Four simulations were combined, each with 1/4 the width of that of llnaba et al.l (|200ll ) at the 
nominal resolution width (Eq. [5|). 



level. Perhaps dynamical friction in our model seems to act a bit stronger than the A^-body results 
suggest, but we do not consider the offset as critical. 

Note the small plateau of the eccentricity curve of the m = 10^^ g bodies near t ~ 10^ yr. As 
discussed above this is due to the transition to the 3D regime. In addition, our expressions have 
that i?df < R vs in the s.d. -regime, whereas iJjf ~ R^s in the d.d. -regime. Therefore, one can say 
that dynamical friction really starts to operate effectively from this point and this explains why 
the curves lie initially (around t ~ 0) much closer together. We remark, finally, that our treatment 
of the s.d. /d.d. -transition regime is rather crude, which explains the erratic behavior of the curves 
at this point (see also our remarks to ward the end of Appendix IB.3p . In order to achieve a better 
match, SIOO and lOhtsuki et al.l (120021 ) continue to empirically modify their analytical expressions, 
but we consider this beyond the scope of this work. 



3.2. Comparison with llnaba et al.l (120011 ) 



Ne xt, we compa r e the outcome of our numerical model including accretion with the statistical 
study of llnaba et al.l (|200ll ) (101 in this section). 101 presents a fully-consistent and highly accurate 
model, includ ing both the dynamica l as we ll as the collisional evolution. 101 in turn compares their 
work against IWetherill and StewartI (j 19931 ) such that we in fact compare three statistical models 
against each other, see Fig. 9 of 101, where the cumulative number of bodies and the velocity 
structure are plotted. 
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The simulation parameters are the following (after Table 1 of IWetherill and StewartI Il993l ) . 
The initial {t = 0) distribution of bodies is monodisperse of mass mo = 4.8 x 10^^ g and horizontal 
velocity v = 4.7 x 10^ cm s~^. The internal density of the bodies is fixed at /)s = 3 g cm~^ and 
the gas density is pg = 1.2 x 10~^ g cm~^, conditions that correspond to a disk radius of 1 AU. 
The distribution is evolved until t = 1.5 x 10^ yr by which time the bodies have grown to masses 
of m ~ 10^^ g, corresponding to sizes of ii 2 X 10^ km. 



101 and IWetherill and StewartI (|l993l ) model a patch of the disk of 0.17 AU. From Eq. [5] we 
obtain Aa = 6.3 x 10~^ AU as the nominal resolution width, which means that 270 zones must be 
included. This proved to be a bit too demanding for the program since each zone must contain a 
minimum number of representative bodies. As a solution, we have reduced the number of zones 
by a factor four {N^o = 67) and computed four of these runs. Each simulation then models a 
total width 1/4 of that of 101, with the four (independent) runs being combined for the total. The 
reduction in simulation width can be justified since the 67 zones is still larger than the minimum 
(Eq. [6]); i.e., the simulation width is still sufficiently wide to harbor R = 2x 10^ km bodies towards 
the end of the simulation. 

The number of particles per zone equals A'res = 600 and the total number of representative 
bodies per run is A'rb = -^rcs-^zo ~ 40000. Other combinations of N^oiNj-^s and Aa will be 
investigated in the next section. 

Figure [6] shows the results. Like Fig. 9 of 101 we show the cumulative number distribution, 
i.e., the number of bodies of mass less than m, and the velocity distribution. For the latter we 
have binned the RBs by mass and shown the mass-weighted planar velocity component. We plot 
the distributions at the same times as in 101. Note that the curves in the velocity plot are in our 
case occasionally a bit noisy, due to the low-N statistics. The gaps in the velocity plot are caused 
by the absence of RBs at these masses. 

Comparing the figures, we find an excellent match. The general trends and shapes of the 
curves are in agreement. There are minor differences but we do not think this invalidates either 
model. Perhaps the biggest difference is the speed of the initial evolution, which by comparing the 
t = 10^ yr curves, can be seen to be faster in our case. During the initial stages growth is rapid (i.e., 
runaway growth) and the high-mass end of the distribution is a rather sensitive function of time, 
see Sect. 14.11 Therefore, we do not consider the offset as critical to the validity of either model. 

We conclude that with the choice of these values for NzojNj-es, and Aa we obtain a satisfactory 
match. Next, we will test how sensitive the results are upon variation of these control parameters. 



3.3. Convergence tests 



Table [3] lists 14 simulation runs where the control parameters N^o, ^res, and Aa are varied. 
The physical parameters like the gas density are kept the same as in the previous section. When 
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Table 3: List of simulation runs to test the influence of the control parameters. 



id 


AT 

A'zo 


AT" 


iVrun Aa 








Irg 






J 2k 


j,min 
<Pvs,50 












AU 






ATT 
AU 


yr 






yr 






(1) 


(2) 


(3) 


(4) 


(5) 






(6) 


{') 






(8) 


(9) 




-1 
1 


4 


zoOO 


1 


4.2 X 


10 


-2 


0.17 


l.z X ID 






0.0 X lU 


o ( \ - y in — 2 

2.9 X lU 




I 


4 


5000 


1 


4.3 X 


10 


-2 


0.17 


1.1 X 1(J 






4.2 X 1(J 


4.7 X lU 




o 

6 


17 


590 


1 


1.0 X 


1 n- 

10 


-2 


0.17 


1 o \ / 1 a3 
l.z X 1(J 






y.D X ID 


1 o \ ^ 1 n— 1 
1.2 X lU 




A 

4 


17 


1180 


1 


1.0 X 


1 n- 

10 


-2 


0.17 


11 \ / 1 n3 
1.1 X ID 






y.u X ID 


ID. (J X ID 




5 


17 


2400 


1 


1.0 X 


10 


-2 


0.17 


9.8 X 10 






1.2 X 10 


11 ^ , 1 r\ — 1 

1.1 X 10 




R 
L) 


U 1 




1 


2.5 X 




-3 


n 1 7 


1 1 V ID'' 








Q n V 1 

O.U A ±u 




7 
i 


1 


300 


1 


2.5 X 


10" 


-3 


0.17 


1.2 X 10^ 






1.4 X 10^ 


3.3 X 10-1 




8 


67 


600 


1 


2.5 X 


10" 


-3 


0.17 


1.2 X 10^ 






1.6 X 10^ 


3.3 X 10-1 




9 


67 


150 


4 


6.3 X 


10 


-4 


0.042 


(1.2 ±0.1) 


X 


10^ 


(2.4 ±0.3) X 10^ 


1.0 ±0.1 




10 


67 


300 


4 


6.3 X 


10 


-4 


0.042 


(1.3 ±0.0) 


X 


10^ 


(2.2 ±0.4) X 10^ 


(9.2 ±0.7) X 


10-1 


11 


67 


600 


4 


6.3 X 


10 


-4 


0.042 


(1.3 ±0.1) 


X 


10^ 


(2.4 ±0.4) X 10^ 


(9.4 ±0.4) X 


10-1 


12 


251 


40 


4 


1.7 X 


10" 


-4 


0.042 


(3.0 ± 1.4) 


X 


10^ 


(2.6 ±0.3) X 10^ 


2.8 ±0.4 




13 


251 


80 


4 


1.7 X 


10" 


-4 


0.042 


(1.2 ±0.0) 


X 


10^ 


(2.4 ±0.2) X 10^ 


2.4 ±0.4 




14 


251 


160 


4 


1.7 X 


10" 


-4 


0.042 


(1.3 ±0.1) 


X 


10^ 


(2.3 ±0.1) X 10^ 


2.4 ±0.2 






Note.- 


— Columns denote: (1) 


simulation identifier; (2 


) number of zones; (3) number of 





representative bodies per zone; (4) number of runs; (5) width of a single zone (6) total simulation 
width; (7) runaway growth timescale; (8) time to grow to 2000 km; (9) minimum filling factor of 
bodies that make up 50% of the viscous stirring power of a zone. The runaway growth timescale 
Trg is defined in Eq. [18l For simulation-id 9-14 where multiple runs have been performed the 



spread in Tj-g and is also given. 

varying the parameters, we have kept the total number of RBs {=N2o x -/Vres) approximately the 
same; that is, an increase of A^^^o by a factor of four is accompanied by a decrease of Nj-es by the 
same factor but for a given N^o we run several models at different A^res- In Table [3] the runs are 
listed by increasing number of zones, Azo- 

Because of the statistical noise associated with the Monte Carlo method, it is rather difficult to 
conduct an unambiguous test for convergence. Rather than focusing on a single measure, perhaps 
the best approach is to compare the evolution of the curves over an extended period. This is done 
in Fig. [7] where the radius of the largest body, Ri, is plotted against time for all 14 runs listed in 
Table [3l The runs are identified by their A^^o with runs having the same A^o sharing the same line 
style. For example, the two A'zo = 4 runs (^^1 and 2) are both identified by a dark grey line. For 
the A^zo = 67 runs we distinguish between two values of Aa, 2.5 x 10-^ and 6.3 x 10"^ AU, the 
lower ones being run at a higher resolution (HR). 

From Fig. [7] we conclude the following: 
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Fig. 7. — Convergence test. The evolution of Ri{t), the radius of the largest body at time t, is 
shown for all 14 simulations of Table [3l Simulations are identified by their number of zones, A'zo- 
The inset shows a zoom in for 1 x 10^ < Ri{t) < 2 x 10^ km. 

1. One of the N^o = 251 simulations {light grey curves) shows very erratic behavior. This 
concerns the Nj-^s = 40 run 12). The low number of RBs per zone turns out to be too few 
to accurately resolve the mass distribution. 

2. Towards the end of the simulation a clear systematic divergence among the curves can be 
seen with the runs of the smallest Aa evolving much faster than those of large Aa. In models 
that do not resolve the disk spatially, the disk appears more homogeneous (in phase space) 
than it actually is. In the low-A'zo models the stirring of planetesimal bodies by the runaway 
body/oligarch occurs at a slower pace because it has to stir so many of them. This artificially 
enhances the accretion. 

3. However, the N^o = 251 (# 12-14) and iV^o = 67 HR (high resolution, # 9-10) curves do not 
separate towards the end of the simulation. This justifies our criterion for Aa, Eq. O 

4. At earlier times these trends are not so obvious. In Fig. [7] the inset shows a zoom of a 
region around Ri ~ 10^ km. The final divergence is not apparent here with stochastic 
behavior seeming to dictate the overall behavior. All curves lie pretty close together, except 
for simulation #12. We also find no clear dependence on A^res- 

To make these findings more quantitative, we have in Table [3] included the timescale to produce 
a 2 X 10^ km body, T2k, and the runaway timescale T^g. The latter is a fit over the exponential 
region of the curve in Fig. [71 i.e., for times t < 10^ yr (see below. Fig. [8)3 and Eq. [T8]) . Except for 
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run 7^12, we do not see a clear trend of Tj-g over the various simulation runs. But for this trend 
becomes obvious. 

In Col. (9) of Table [3] we have listed the minimum filling factor (/"^lo massive bodies 

that together make up more than 50% of the stirring power in a certain zone. This quantity is 
obtained as follows: 



1. We select the masses above mvs,5o which together comprise 50% of the viscous stirring power. 
Since the stirring power of the bodies scales as Sm (Eq. [T]), i.e., as the squares of the masses 
of the bodies, the criterion for mvs,50 is 

mf/j;m2 = 0.5. (15) 

mi>mvs,50 / i 

2. We sum up the 'spheres of influence' concerning accretion and compare this to the width of 
the zone Aa. Here, we assume that the bodies are dynamically cold and that the net impact 
radius is on the order of the Hill sphere. Thus, 

(t>vs,50 = yZ '2Rh{mi), (16) 

mi>mvs,50 

3. We take the minimum of </>vs,50 over the course of the simulation run to arrive at i?5>™™o 



During the runaway growth phase the massive bodies become dynamically very cold, moving 
on circular orbits. When scattering dominates we may expect the bodi es to become isolated s ince 
scattering, together with dynamical friction, leads to orbital repulsion (jKokubo and Idalll995l ). If 
the stirring power resides in the (few) big bodies, their mutual scattering may then cause them 
to become isolated in the sense that their mutual spacing in terms of semi-major axis becomes 
too large for collisional interactions. Then, i;^vs,50 drops below unity, in which case our statistical 
assumption ~ that the bodies are uniformly distributed over the width of the zone - breaks down. 
However, when i;^vs,50 ^ 1 scattering among the bodies cannot result in their isolation; the viscous 
stirring power is shared among a sufficiently large number of bodies to warrant the validity of the 
statistical assumption. 

We find that, initially, 4'vs,50 S> 1 since the stirring power is initially determined by the plan- 
etesimals (mvs.so ~ "^o)- However, during the fast runaway growth phase a power-law distribution 
emerges, in which the stirring power becomes dominated by the bodies at the high-mass end; (/'vs,50 
then quickly decreases and a minimum is obtained. In Table [3] the minimum of 0vs,5O is given. We 
see that for the large Aa runs it falls below unity, indicating that the simulations do not resolve 
the spatial structure properly. Indeed, as can be seen from Fig. [71 growth in the low N^o is the 
fastest, but this growth - triggered by merging of big bodies - is artificial. However, in the large 
Nzo models 0vs,5O stays above unity, indicating that the simulation is properly resolved. 
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Fig. 8. — (left) The surface density spectrum, dT,/dR for the 1 AU, gas drag simulations discussed 
in the previous section. Curves are plotted at every factor of 2 increase in the radius of the biggest 
particle, Ri{t). The dotted line in the upper- left of the panels shows the trend of the column density 
spectra Ns(m) oc with p = —2.5. The dashed line in the bottom-right corner corresponds to 1 
particle per bin. {right) Several statistics shown as function of the evolutionary parameter Ri{t): 
(i) the ratio of Mi to the characteristic mass of the population, m*; (ii) the maximum velocity 
within the population to the Hill velocity of the largest body; and (iii) the time (on the second 
y-axis). The dotted vertical line denotes the minimum of Vx/vh- 

These findings indicate that the criteria outlined in Sect. [231 i.e., Eqs. [5] and El do lead to 
numerical convergence. For the number of RBs per zone we recommend at least A'res = 100. 
However, we will typically use a larger A'res to reduce the MC-noise. 



In Fig. [8] we again plot results from the previous section - i.e., run #11 of Table [3l which 
results were also presented in Fig. [6l Figure [8^ now shows the mass spectrum, instead of the 
cumulative distribution of Fig. [6^. (Thus, {dT, / dR) A.R gives the surface density of bodies within 
the size interval [R, R + AR]). For our purposes we find it more instructive to show mass spectra 
like Fig. [8^ rather than cumulative distributions, since the relevant features turn out more clearly. 
However, Fig. [6ti is in fact just an integrated copy of Fig. [8^. 

Once the radius of the maximum particle in the distribution has increased by a factor of 2, 
a curve is plotted and the corresponding time is indicated. In Fig. [8^ the dotted auxiliary line in 
the upper left corner indicates the trend if the column density spectrum, Ns{m) = {\/m)dTi/ dm, 
would be a power-law of mass, Ns(m) oc mP with exponents p = —2.5. For reference, a flat slope 



4. Runaw^ay growth vs. oligarchy 



4.1. Runaway growth indicators 
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{p = —2) would indicate that the distribution contains an equal amounts of mass per logarithmic 
size bin. Thus, Fig. [8^ shows that the high mass tail of the distribution first flattens towards a 
p ~ —2.5 slope, before it breaks, and that most of the mass remains at the initial planetesimal 



size R o. This behavior is consistent with A^-body simulations (|Kokubo and Idalll996l : iBarnes et al 



20091 ). Afterwards, the distribution evolves into two components with two identifiable peaks: one 
bump appears at low-masses that with time evolves to larger sizes and another 'spike' is associated 
with the largest bodies. Finally, the line in the lower right corner of Fig. [8^ indicates the size 
spectrum if there would be only a single body in the mass bin. The curves have to stay above this 
line. 

What would be the best indicator to assess whether a system is in runaway growth (RG)? 
Probably the best indicator is AI1/M2, i.e., the ratio between the mass of the biggest to the second- 
biggest body in the system. When this ratio increases, the system is in RG; otherwise it is not. 
Unfortunately, the problem is that this quantity behaves very erratically when the bodies are still 
close to each other in terms of their masses: stochastic processes then interfere to produce a noisy 
behavior. For this reason, rather than M1/M2, we propose to use the ratio of the most massive 
body over the 'characteristic mass', Mi/m*, as an indicator for RG, where m^: is defined as 

= ^ N{mi)mf / ^N{mi)mi , (17) 

which traces the particles that contain most of the mass, excluding the most massive body Mi. For 
narrow distributions m* approximately corresponds to the peak of the m^Nirn) mass spectrum. 
However, here we will merely use it as a tracer for the 'mass flow'. If increases more steeply 
than Ml, the mass flow is no longer preferentially directed to Mi; the accretion rate of other bodies, 
e.g., m ~ M2 or m ~ m*, then starts to outweigh that of the most massive one. In two component 
systems, pr evious toy models have shown that Mi /m^, increases for RG; o therwise the system is 



not in RG (Wetherill and Stewart 1989; Lee 



200c 



Ormel and Spaanal2008l ). 



This ratio is plotted in Fig. [S}3 by the solid black line. There are still many stochastic fluc- 
tuations due to merging of large bodies, particularly at later times, but the general trend is clear. 
At Ri ~ 100 km Mi/m* reaches a plateau and starts to decline more visibly after Ri ~ 400 km, 
indicating that RG has terminated. The mass of the system flows to larger sizes, but not exclusively 
to one object. 

This behavior can be understood from the trend of Vx/vh^ plotted by the dashed curve. Here, 
Vx is the maximum random velocity in the simulation (which is associated to the low-m planetesimal 
bodies from which the runaway body is accreting) and Vj^^ is the Hill velocity of the most massive 
particle. Thus, the ratio Vx/vh is a measure of the amount of gravitational focusing (GF) the 
runaway body experiences when accreting the small bodies: if it decreases, GF increases, whereas 
if Vx/vh increases, GF decreases (recall from Sect. 12. ll that relates to the escape velocity as 
Uesc = Vh^/oja). If only accretion would operate (constant Vx) Vx/vh oc R^^ . However, by exciting 
the random motions of the planetesimals, viscous stirring counteracts the decrease of Vx/v^: it 
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increases Vx and decreases the GF. Prom Fig. [8)o it can be seen that accretion dominates during 
the initial stages but also that the decrease of Vx/vh is not so steep as in the stirring- free limit. 
We find that the scaling is now rather Vx/vh oc R^^'^ . Clearly, there is a positive feedback effect 
at work: the growth of the largest body increases Vh, which in turn increases the GF factor. This 
mechanism operates in the initial stages. However, the relative importance of the viscous stirring 
becomes more apparent at low Vx/vh-, see Fig. [31 Consequently, a minimum of Vx/vh is reached. 
We denote this point {Ri ~ 400 km) the transition size i2tr- GF factors peak at Ri = i?tr- 

The transition signifies a different growth phase as can be seen by the gray curve in Fig. [8)3, 
which shows the simulation time as function of Ri. At times Ri{t) < i?tr growth proceeds expo- 
nentially, Ri{t) oc exp[t], whereas if Ri{t) > i?tr the growth proceeds much slower. Note the linear 
spacing of the time-axis in Fig. [Sb- We obtain the associated timescale T^g empirically by a fit to 
Ri{t), see Fig. [Hb, i.e., 

Mi(t)«exp(^^) . (18) 

(Note that T^g is defined in terms of mass, not radius). We refer to Tj-g as the runaway-growth 
timescale since it is associated to the RG part of the evolution; i.e., the phase where R < -Rtr- 

The observation that the growth before the transition point proceeds exponentially at a 
'measured' rate Tj-g is, above all, an empirical finding. It implies that the accretion rate of 
the largest particle dMi/dt proceeds linearly with its mass Mi. From the observed relation 
Vx/vh OC Ri ^^"^ oc one indeed finds that the accretion rate in the dispersion-dominated 

regime is (approximately) linear, dMi/dt oc R\{vh/vx)'^ oc Mi. In Sect. I5.3T] we will see that dur- 
ing the runaway growth phase all size bins contribute (approximately equally) to the growth of the 
biggest body. Mergers among big bodies take place in the s.d. -regime, for which the accretion rate 
is also linear with the mass of the biggest body. 



4.2. The emergence of oligarchy 

Next, we consider the evolution beyond the runaway growth transition Ri{t) > i?tr- Figure [9] 
shows the spatial distribution of the planetesimal swarms corresponding to the simulation discussed 
above, i.e., model #11 of Table [3l Each symbol indicates a planetesimal swarm (representative 
body). In Fig. [9] the total mass in the planetesimal swarms is indicated by the area of the symbols. 
Representative bodies that contain a single particle {Ng = 1) are indicated by diamonds rather 
than dots. The Hill spheres of the ten most massive bodies are indicated by red bars. Finally, the 
colorbar represents the gravitational focusing factors like in Fig. [SJ the ratio of the random velocity 
(v) to the Hill velocity of the largest body (vh)- 

Figure [9] shows the planetesimal distribution at three different times: at Ri = 500 km, just 
after the transition size; at -Ri = 10^ km; and at i?i = 2 x 10^ km, the end state of the simulation. 
Figure [9] very clearly shows that with time: (i) the gravitational focusing factors increase; (ii) the 
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Fig. 9. — Scatter plot of bodies' mass and position at three times during the ohgarchy phase, 
corresponding to the lAU gas-drag simulation. Each dot represents a planetesimal swarm. The 

1/3 

size of the dot is a measure for the total mass of the swarm, such that area(dot) oc ms4arm- 
Individual bodies (i.e., those that have Ng = 1) are shown by diamonds. The color is a measure for 
the random velocity (eccentricity) of the bodies in terms of v/vh (see the colorbar on the right). 
These values are normalized to the Hill velocity (vh) of the largest body, which value is indicated 
in the panel. The Hill radius of the 10 most massive bodies are indicated by a red bar. Note that 
the scaling on the x-axis differs among the three panels. 



number of oligarchs decreases; (iii) the gap in mass between the oligarchs and (leftover) planetesimal 
increases. The oligarchy gets more pronounced with time. 

For Ri{t) > Rtr, viscous stirring by the high-mass bodies gains the upper hand and accretion 
times increase due to the fact that gravitational focusing factors decline (increasing Vx/vh)- Other 
(massive) bodies then catch up. Indeed, Rtr also approximately corresponds to the point where 
Ml / m^, starts to decrease noticeably, see Fig. [51 In due time the gas drag should balance the 
stirring to produce an equilibrium eccentr icity that is characterized by a constant v^/vh (e.g.. 



Ida and Makindll993l : lThommes et al.ll2003l ). Our simulation evolves towards this state in a rather 
erratic way that is caused by the merging among big bodies. Note that in our case Vx traces the 
largest random velocity of any body. Towards the end of the simulation, in Fig. [9]3, a reversal in the 
random velocity distribution has taken place, in which the low-mass planetesimals have no longer 
the largest random velocity due to the fact that gas drag is more effective for these bodies (see also 
Fig. E]). This explains why Vx/vh does not readily approach a constant value. 

Since the number of oligarchs declines, a large contribution to the growth during the oligarchy 
stage should come by the merging of similar-size bodies. The oligarchs are simply packed too 
densely to guarantee their mutual existence. There is some limited diffusion of the oligarchs due to 
merging of bodies (see Sect. [2r3]) . Scattering among oligarchs is not implemented in our approach. 
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However, even if implemented, scattering cannot isolate these bodies from each other since there 
are too many of them. The only way forward is to merge the oligarchs. The consequence is that the 
distance among the oligarchs in our simulations is therefore at least several Hill radii - a spacing 
that directly reflects the choice of Rco\ = 2.5Rh in the low velocity regime. Oligarchs that come 
within this distance have a strong probability to merge. Although we cannot reproduce the exact 
spatial structure of the A^-body simulations, this pictur e of merging oligarchs is broadly consistent 
with the A-body simulations of iKokubo and Idal (|2002l ). 



Since our model does not treat scattering, gap formation within the planetesimal disk is not 
observed. This effect coul d serious l y slow down the growth during the oligarchy stage. Judg- 
ing from the calculations of iRafikovl (|2003l ). gap formation becomes important for bodies of mass 
> 10^^ g (i? « 10^ km). Yet, again, the oligarchs are likely to be too densely packed to prevent 
them accreting planetesimals. With increasing size of the oligarchs scatterings should become more 
pron ounced, however. In a recent A'^-body simulation involving Earth mass embryos and planetesi- 
mals, iLevison et al.l ()201Cl ) observed that the planetesimals were scattered over AU-distances - out 
of the feeding zones of the embryos. Modeling these effects are beyond the scope of this work. 



4.3. Definitions of runaway growth: local and global 

A system of bodies is in RG when the ratio M1/M2 increases with time, i.e., d{AIi / M2) / dt > 



(jWetherill and Stewar till 9891 ). where Mi and M2 are, respectively, the mass of the most massive 
and the second-most massive particle in the system. Alternatively, we can compare the accretion 
timescales; thus, the system is in RG at time t when the condition 

Tf"(t) < T|"(t) (19) 

is fulfilled, where T'^'^ is the accretion timescale 



1 dM^ 



M dt 



(20) 



In a two-component system, the mass accretion rate dM/dt is proportional to the collisional 
cross section i?^oii approach velocity Va, and the number density of particles Ns/2hes, i.e., 
dM/dt cx R'^^iVaNs/heg. If we take the dispersion-dominated regime where i?coi R'^{vh/va)'^ , 
Va/hes = ^ (see Sect. 12. 4j ). and take Nsi = Ns2, the condition Eq. [T2] translates into 

i?,f^y<i?jMy. (21) 

\VhlJ \Vh2J 

If the interactions take place in the same zone, i.e., Vxi = Vx2, the RG-condition is always satisfied 
since Vh oc R. In that case the RG-index k as in 

^ocM^ (22) 
dt ^ ' 
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equa ls k = 4/3. The usual criterion is that k > 1 is needed for RG to ensue (jWetherill and Stewart 



19891 ) ■ However, this is only valid for interaction within the same spatial zone, i.e., where competitor 



bodies accrete from the same reservoir of (low-mass) planetesimals. 

Similarly, if Vx is constant everywhere the RG-condition holds globally. However, this is not to 
be expected since viscous stirring will increase the random velocities of the planetesimal bodies. We 
saw above that initially, growth proceeded exponentially, Vx/vh oc R\ with q ~ —0.5, which means 
that the RG-condition is marginally satisfied. Nonetheless, we will refer to the initial (exponential) 
growth state as the runaway growth phase. After the transition size has been reached {Ri > i?tr) 
Vx/vh increases with Ri and Eq. [19] is no longer globally fulfilled (but locally it still is since Vxi = 
Vx2)- This, is the oligarchy stage. 

As the term 'runaway growth' is sometimes used rather colloquially in the literature, we sum- 
marize a few statements that are in agreement with the formal definition above, Eq. [T9l 

1. RG is not always synonymous with fast growth. Runaway growth is, for understandable 
reasons, often identified with fast growth rates, and, therefore, sometimes associated with the 
s.d. -regime, Vx < Vh- However, in the s.d. -regime k < 1 and the growth mode is not runawa3ll|. 
Indeed, Eq. [19] does not specify the absolute growth rate of the bodies. 

2. RG does not require dynamical friction. During RG bodies interact in the d.d. -regime. When 
bodies start out at v > v^sc no RG takes place. In this case dynamical friction, which reduces 
the random velocities of the most massive bodies, could shift the most massive bodies in the 



d.d. -regime, such that Eq. [19] materializes (jWetherill and Steward! 19931 ). However, dynamical 



friction is not required to sustain RG as long as they take place in the d.d. -regime. 

RG is not synonymous with immediate mass separation. Although Eq. [T9] implies that masses 
will separate, this criterion neglects stochastic effects and the actual mass-doubling time of 
the second-most massive body can still be shorter. In the long run, however, RG will result 
in mass-separation but Eq. [19] does not specify at which point this occurs. 

In oligarchy, the system is only locally in RG. At late times, the runaway body (or oligarch) 
regulates the velocity dispersion of the planetesimals leading to a positive feedback on the 
random velocities of the bodies, which slows down the growth. As explained above, inter- 
actions within the same spatial zone, i.e., the region of the disks where the runaway body 
have stirred the random velocity of the planetesimals to (the same) Vx, always fulfill the RG- 
condition as long as the d.d. -regime holds. However, among bodies of different spatial zones 
Eq. [19] is no longer satisfied. The c ombined effect of local RG and isolation is better known 



as oligarchy (jKokubo and Idalll998l ). 



^For a thick planetesimal disk in the shear-dominated regime we find dM/dt oc A4 and k = 1; for a very thin 
planetesimal disk, on the ot her hand, flcoi becomes larger than the scaleheight /left. In that essentially 2D setting 
dM/dt oc M^/^ and /t = 2/3 (jKokubo and Idall 19961 ). 



-37- 



In the remainder of the paper we will identify systems that experience the initial exponential 
(global runaway) growth with the runaway growth phase and systems that undergo only local-RG 
with the oligarchic growth phase. That is, we consider an evolutionary sequence in which the 
oligarchy phase supersedes the runaway growth phase. We distinguish the runaway growth and the 
oligarchy phases of the planetesimal accretion process as follows (see Fig. [Sb): 

• In the RG phase {Ri{t) < Rtr), Vx/vh decreases and GF- factors increase. Growth occurs 
exponentially at a characteristic timescale T^g. In addition, the quantity Mxjm^ increases 
during most of the phase. 

• In the oligarchy phase {R\{t) > Rtr), Vx/v^ increases (GF-factors decrease). As a result, 
accretion timescales increase rapidly and Mi/m* decreases. 

5. Parameter study 

Table H] contains the list of runs that have been performed. The prime goal is to cover a wide 
range of physical conditions and disk radii to see whether the picture sketched in the previous 
section holds generally. In particular, this concerns the behavior of the Vx/vh indicator and its 
relation to the growth rate and the Mi/m* indicator. For this reason we perform simulations at 
three disk radii: 1, 6, and 35 AU, in which damping by gas drag and fragmentation are varied (first 9 
entries in Tabled]). Since we find that fragmentation is an important mechanism, we next test how 
sensitive the outcome is under variation of the fragment size Ofi- and the coefficient of restitution 
parameter, e (see Sect. [23]) . The next class of runs contain turbulent stirring, abbreviated Ts. 
Finally, we focus on the gas-free runs at 35 AU and vary additional physical parameters, like the 
initial planetesimal size -Rq, the initial random velocity fo, and the surface density S. 

In Table H] Col. 1 gives the model name, which is a mnemonic abbreviation of the semi-major 
axis at which the run was performed and the features it includes. Column 2 lists the number of 
superparticles (swarms) used per zone and Col. 4 the number of zones that are included, which 
follows the guidelines outlined in Sect. 12.31 Column 4 provides the surface density of solids. Col. 5 
lists the total width of the simulation. Col. 6 the initial planetesimal radius, and Col. 7 the initial 
random velocity dispersion, which is taken to be half the initial escape velocity of the bodies. Bodies 
corresponding to simulations at 1 AU have an internal density {ps) of 3 g cm~^, whereas the bodies 
from the 6 and 35 AU simulations have an internal density of 1 g cm~'^. 

The model feature abbreviations (Col. 1) imply that simulations include: 

• Velocity evolution (Ve). This is shorthand for all processes through which the random veloc- 
ities are affected by collisions and gravitational stirring: viscous stirring, dynamical friction, 
and collisional cooling. Inclinations are calculated independently from eccentricities. All runs 
include these features. 
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Table 4. List of simulation runs 



Radius+Features 


^ * res 




E 


Aa 


Rq 


vo 


Comments/Fij 


jure refs 


[AU] 






[g cm"2] 


[AU] 


[km] 


[m/s] 






(1) 


(2) 


(3) 


(4) 


(5) 


(6) 


(7) 


(8) 




Models at different disk r 


adii 














fVe 


500 


67 


16.8 


0.042 


7.3 


4.7 


Fig. [10] 




fVeGd 


500 


67 


16.8 


0.042 


7.3 


4.7 


Fig. [10] 




IVeGdFr 


500 


67 


16.8 


0.042 


7.3 


4.7 


Fig. [TO] 




6Ve 


500 


33 


2.0 


0.18 


7.3 


2.7 


Fig. [TO] 




6VeGd 


500 


33 


2.0 


0.18 


7.3 


2.7 


Fig. [TO] 




eVeGdPr 


500 


33 


2.0 


0.18 


7.3 


2.7 


Fig. [TO] 




35Ve 


500 


13 


0.20 


0.97 


7.3 


2.7 


Fig. [TO] 




35VeGd 


500 


13 


0.20 


0.97 


7.3 


2.7 


Fig. [TO] 




35VeFr 


500 


13 


0.20 


0.97 


7.3 


2.7 


Sect. [O] Fig. 


Ho] 


Models varying fragmentation parameters 












35VeFrHi-afr 


500 


13 


0.20 


0.97 


7.3 


2.7 


flfrag = 10 cm 




35VeFrHi-e 


500 


13 


0.20 


0.97 


7.3 


2.7 


e = 0.1 




35VeFrLo-e 


500 


13 


0.20 


u.y ( 


i .O 


Z. i 


e — U.UUi 




Models including 


turbulent stirring 














IVeGdFrTs 


500 


67 


16.8 


0.042 


7.3 


4.7 


Fig.[n] 




6VeGdFrTs 


500 


33 


2.0 


0.18 


7.3 


2.7 


Fig.[n] 




35VeGdFrTs 


500 


13 


0.20 


0.97 


7.3 


2.7 


Fig.[n] 




35 AU, miscellaneous 
















35VeFrlkm 


500 


95 


0.20 


0.99 


1.0 


0.37 


Fig. [TO] 




35VeFr50km 


500 


5 


0.20 


2.6 


50.0 


19.0 


Fig. [TO] 




35VeFrLo-vO 


500 


13 


0.20 


0.97 


7.3 


1.0 


Fig. [TO] 




35VeFrHi-vO 


500 


13 


0.20 


0.97 


7.3 


25.0 


Fig.[lS] 




35VeFrHi-I] 


500 


13 


2.0 


0.97 


7.3 


2.7 


Fig. [TO] 




35VeFrLo-E 


500 


13 


0.020 


0.97 


7.3 


2.7 


Fig. [TO] 





Featured abbreviations, which make up the simulation name listed in Column 2, denote: Ve, 
velocity evolution (includes viscous stirring, dynamical friction, and coUisional cooling); Gd, gas 
drag; Fr, fragmentation; Ts, turbulent stirring, Lo-vO, low initial random velocity, Hi-vO, high initial 
velocity, etc. 

Note. — Columns denote: (1) disk position (in AU) and features; (2) number of simulation 
particles per zone; (3) number of zones; (4) surface density in solids; (5) total simulation width; (6) 
initial radius of bodies; (7) initial random velocity; (8) comments. 
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• Gas drag (Gd), which damps the random velocities according to Eq. [TTJ 

• Fragmentation (Pr), for approach velocities that exceed Vesc (i-e., the case without GF), 
according to the procedure outlined in Sect. 12. 5[ 

• Turbulent stirring (Ts), an external source of excitation to the random motions of the bodies, 
see Sect. 12.5.31 We list only runs with 7 = 10~^. Runs that were performed at 7 = 10~^ did 
not result in accretion. 

Simulations are continued until a size R = 2 x 10^ km is reached for the most massive bodies. 
Table [5] provides statistical quantities that characterizes the outcome of the runs. These include, 
the radii at which Mi/m* and Vx/vh have their extrema, denoted i?* and Rtr, respectively, and 
the runaway growth timescale, T^g (see Sect. I4T]) . When error bars are given these indicate the 
spread in the quantities over 5 independent simulation runs. Table [5] further gives the mass fraction 
fragments constitute at the end of the simulation, and the total fraction of the mass that has 

By comparing these fractions one gets an estimate of the importance 
of fragmentation and of the importance of the re-accretion of these fragments. 

5.1. Runs including gas drag and fragmentation 

Figure [TOl presents the quantities Mi/m*, Vx/vh, and time for a variety of simulations including 
gas drag and/or fragmentation. Runs at three different disk radii are shown: 1 AU, 6 AU, and 35 
AU. Note that in Fig. [TOb time is plotted in units that depend on the radial distance: Myr for the 
35 AU runs, 10*^ yr for the 6 AU runs, and 10^ yr for the 1 AU runs. 

The general trend of the curves is the same for all the disk radii and reflects the discussion 
of Sect, im First, in the runaway growth phase, Mi/m^, rises and Vx/v^ decreases. The outer 
disk models start out at a larger value of Vx/vh since Hill velocities for the same mass decrease 
with increasing disk radius. Very generally, one can say that interactions in the inner disk are 
more prone to take place in the shear-dominated regime. At the point where Vx/vh has reached a 
minimum, accretion timescales increase. A major difference is that for the 1 AU run the transition 
occurs much sooner (in terms of the evolutionary parameter Ri{t)) than for the outer disk models. 

From Fig. [10] it is seen that gas drag does not influence the growth during the initial runaway 
growth phase (i?i < -Rtr), indicating that it acts on a longer timescale. The exception is the 35 AU 
model, where the divergence occurs at a relatively low Ri. The likely reason for this phenomenon 
is the different gas drag law {Co > 1) small bodies experience at 35 AU, which, somewhat para- 
doxically, increases the (relative) effectiveness of gas drag. However, in our simulations we do not 
reduce the gas density at large timescales; for t ~ 10^ yr the gas would surely have dissipated from 



Formally, this number can exceed unity since we do not adjust for multiple fragmentation-accretion cycles. 



Table 5. Simulation results 



Name 



(1) 



[km] 
(2) 



-Rtr 

km] 
3) 



[Myr] 
(4) 



T2k 

[Myr] 
(5) 



fend 
./frag 

(6) 



fend 
^tot 



(7) 



IVe 

IVeGd 

IVeGdPr 

6Ve 

6VeGd 

eVeGdPr 

35Vc 

35VeGd 

35VePr 

35VePrHi-afr 

35VeFrHi-e 

35VeFrLo-e 

IVeGdPrTs 

6VeGdFrTs 

35VeGdFrTs 

35VeFrlkm 

35VeFr50km 

35VcFrLo-vO 

35VeFrHi-vO 

35VeFrHi-E 

35VeFrLo-S 



(1.9 ±0.4 
(2.0 ±0.3 
(1.9 ±0.2 
(3.4 ±0.3 
(3.5 ±0.8 
(3.5 ±0.5 
(5.0 ±0.5 
(5.4 ± 1.1 
(4.5 ±0.3 
(4.7 ±0.9 
(2.7 ±0.8 
(4.9 ±0.7 
3.5 X 10^ 
7.1 X 10^ 

1.1 X 10^ 
1.7 X 10^ 

1.2 X 10^ 
5.0 X 102 
3.4 X 10^ 
7.9 X 10^ 
2.4 X 10^ 



X 10^ 
X 10^ 
X 10^ 
X 10^ 
X 10^ 
X 10^ 

X 102 
X 102 
X 102 

X 102 
X 102 
X 102 



3.8 ±0.7 

3.6 ±0.5 
3.5 ±0.3 
7.2 ± 1.3 

7.7 ± 1.3 
7.2 ±0.9 
1.1 ±0.3 
9.5 ±2.1 
1.4 ±0.3 
1.1 ±0.6 

9.8 ±3.4 

9.9 ± 1.1 
6.2 X 102 

2.0 X 10^ 

2.1 X 10^ 
6.8 X 102 
1.8 X 10^ 

1.7 X 10^ 

1.5 X 10^ 

1.6 X 10^ 

4.8 X 102 



X 102 
X 102 
X 102 
X 102 
X 102 
X 102 

X 10^ 

X 102 

X 10^ 

X 10^ 

X 102 
X 102 



(1.4 ±0.1 
(1.3 ±0.1 
(1.2 ±0.1 
(6.9 ±0.3 
(6.4 ± 0.4 
(5.9 ±0.3 
(1.2 ±0.1 
(5.9 ±0.3 
(5.6 ±0.2 
(1.1 ±0.1 
(1.3±0.1 
(1.2±0.1 
1.0 X 10-2 
9.0 X 10-1 
1.2 X IQi 

6.8 X 10-1 

5.7 X IQi 

5.9 X 10" 
4.2 X 10-1 
4.0 X 10-1 

8.8 X IQi 



10-3 

10-3 
10-3 

10-2 
10-2 
10-2 

lOi 
10° 
10° 
IQi 
10" 
IQi 



(4.1 ±0.4 
(2.3 ±0.2 
(4.8 ±0.2 
(3.0 ±0.7 
(2.1 ±0.1 
(1.2 ±0.1 
(3.2 ±0.7 
(1.6 ±0.1 
(1.1 ±0.0 
(1.7±0.1 
(3.7±0.1 
(2.1 ±0.1 
2.5 X 10-1 
2.2 X lOi 
7.2 X 10^ 

1.8 X lOi 

5.4 X 102 
1.0 X 102 
3.2 X lOi 

8.9 X 10° 

1.5 X 10^ 



X 10-1 
X 10-1 

X 10-2 

X 10° 
X 10° 
X 10° 

X 102 
X 102 
X 102 
X 102 

X lOi 

X 102 













(10.0 ± 0.9) X 10-2 (3.1 ± 0.3) X 10-1 





(4.6 ± 0.9) X 10-2 (6.9 ± 1.3) x 10-2 






(9.3 ±2.6) X 10-3 
(2.4 ±0.4) X 10-2 
(2.2 ±0.4) X 10-2 
(3.4 ±0.3) X 10-3 
3.4 X 10-2 

2.8 X 10-2 

4.4 X 10-2 

3.5 X 10-2 
2.1 X 10-3 
6.7 X 10-3 

1.9 X 10-2 

3.5 X 10-3 

1.6 X 10-2 






(1.7 ±0.4) X 10 
(3.3 ±0.5) X 10 
(3.8 ±0.7) X 10 
(6.7 ±0.6) X 10 
2.0 X 10-1 

2.5 X 10-1 

3.3 X 10-1 
4.8 X 10-2 

4.6 X 10-3 

1.4 X 10-2 
3.0 X 10-2 

6.5 X 10-3 
7.4 X 10-2 



2 

2 
-2 
-3 



Note. — Column entries denote: (1) model name; (2) size at wliicli Mi/m* readies its maximum; (3) size at which Vx/vu reaches its 
minimum; (4) runaway growth timescale; (5) time at the end of the simulation run (at i?i = 2 x 103 km); (6) mass fraction in fragments 
at end of simulation; (6) total mass fraction that has been processed through fragments. Error bars denote the spread over 5 runs. 
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Ri [km] 

Fig. 10. — The statistical quantities for runs at several disk radii. As function of evolutionary 
parameter Ri(t) - the radius of the most massive particle - are shown: (A) Mi/m*, (B) Vx/vh, 
and (C) time. Black curves correspond to runs that include neither gas drag nor fragmentation, 
light gray curves indicate runs that include gas drag, and dark grey curves indicate fragmentation. 
In (C) times are normalized to 10^ yr (for the 1 AU runs), 10^ yr (6 AU), and 10^ yr (35 AU), 
respectively. 



- 42 - 



the nebula. On the other hand, for the longer timescales that characterize the oligarchic stage, 
gas drag increases in importa nce. Balancing stirrin g with damping, it is found that GF-factors 



approach an equilibrium (e.g.. iKokubo and Idall2002l ). Thus, the black curves (no gas drag) should 



separate from the other curves in the oligarchy stage in Fig. llOb . Although we can see the start of 
this process, its overall signature is not very clear yet due to the strong fluctuations in the curves 
caused by merging among similar-size big bodies (cf. also our discussion in Sect. 14.21 

Fragmentation, although not playing a major role in the initial runaway growth phase, has 
the tendency to smooth the features associated with the transition size. There is still a clear peak 
in Mi/m^ but the signatures of the transition are not so evident (but still present) in v^/vh and 
the timescales plot. Compare, for example, the behavior of the dashed curves in Fig. [TOl The 1 
AU and 6 AU curves visibly steepen after the RG/oligarchy transition (at Ri ~ 10^ km) but for 
the 6VeGdFr curve this effect is much less obvious. Fragments do not conform to the self-regulated 
aspect of oligarchy, since their random velocity is not strongly affected by (the growth of) the 
biggest bodies. 

Our simple fragmentation model contains two free (uncertain) parameters that affect the be- 
havior of fragments: the size of the fragments Ofi- = 1 mm and the coefficient of restitution, e = 0.01. 
We have tested the influence of these canonical values by running models with afr = 10 cm and 
varying e with respect to the (gas-free) 35VeFr run, see Table HI A larger fragment size results in 
less efficient cooling among the fragments due to their reduced total cross section for interaction. 
As a result, the fragments are not so efficiently re- accreted during the runaway growth phase and 
the runaway growth timescale, T^g, is rather like the non- fragmentation 35Ve run. The influence 
of varying e is quite significant. This can be understood by considering the extreme limits: e = 
implies that each collision is fully inelastic and no fragmentation takes place, whereas e = 1 implies 
that every v > Vesc collision completely shatters both bodies. Thus, our e = 0.001 run lies closer to 
the non-fragmentation run, whereas in the e = 0.1 run more fragments are produced, resulting in 
shorter overall accretion timescales. 



5.2. The effects of turbulent stirring 

Figure [TT] presents the key indicators for models including turbulent stirring (Ts) and frag- 
mentation (Fr). Models without stirring are plotted for comparison (black lines). Here, only runs 
with a value of 7 = 10~^ for the turbulent stirring parameter are shown since it was found that 
when the stirring parameter 7 was set to 7 = 10"'^ no accretion takes place. In Fig. [TTb times are 
normalized by the fiducial 

and plotted on a logarithmic y-axis. In Eq. [23} tmn is, upon neglect of a numerical constant, equal 
to the initial collision timescale between the bodies of size Rq, internal density ps, and surface 
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^1 [km] 

Fig. 11. — Statistics of runs including turbulent stirring (Ts). The non-Ts runs are plotted for 
comparison (black lines). 

density S. That is, trun = {noaoVQ)~^ with no = T,o/moheS,o, m-o ~ PsRq, ho ~ vo/^, etc. For 
exponential growth or runaway growth one can expect the initial timescale to be a characteristic 
timescale of the system. 

We find that turbulent stirring at fixed 7 affects the outer disk more than the inner disk. 
This can be understood since gas damping is more effective in the inner nebula. As can be seen 
from Fig. [TTl at 1 AU the effects of turbulent stirring are relatively minor. However, adding even 
a relatively small amount of stirring does increase accretion timescales by a factor of 10 or more 
due to the lower focusing factors (Fig. Illb). In addition, the accretion of the biggest bodies is 
dominated by fragments. This has the tendency of erasing the imprints of the runaway /oligarchy 
transition as present in the conventional models (without Ts). This effect was already seen in the 
non-turbulent models, but becomes now more pronounced. 
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Fig. 12. — The evolution of the quantities /frag (the fraction of the total mass that resides in frag- 
ments at that time; solid curve) and /tot (the total mass fraction that has once been in fragments; 
dashed curves). The difference between /tot and /frag is due to accretion of fragments. 

To further illustrate this behavior, Fig. [12] plots the mass fraction in solids as function of 
time, or rather, the evolutionary parameter, Ri{t). Here, /frag is the mass fraction that resides 
in fragments at Ri{t) and /tot the cumulative mass fraction, i.e., the total amount of mass that 
has once been in fragments. The difference between these curves then provides a measure for the 
amount of fragments that are accreted. From Fig. [12] it can be seen that in the turbulent stirring 
run (IVeGdFrTs) collisions quickly put ~1% of the mass into fragments and that these fragments 
are accreted efficiently since /tot keeps increasing while /frag levels out. On the other hand, in the 
IVeGdFr run the importance of fragmentation (and subsequent accretion) only gradually increases. 
It is rather unimportant during the runaway-growth phase, but gathers pace towards the end of 
the simulation, even overtaking the turbulent stirring curves. Therefore, the tabulated fractions 
towards the end of the simulation run (i.e.. Cols. 4, 5 of Table [5|) do not reflect the importance of 
the fragmentation during the earlier runaway growth phase. 



The Ts-runs conducted at 6 and 35 AU (Fig. [TT]) show very interesting behavior. Relative 
velocities are initially so large that GF is unimportant. However, these collisions are just not yet 
violent enough to be fully destructive: there is net accretion and the mass spectrum is characterized 
by an exponential tail at the high mass end. Accretion timescales become very long, though, when 
GF is absent. A self-similar size distribution emerges in which Mi/m* ~ 10. However, at the point 
where vm ~ fesc,M collisions enter the d.d. -regime. This occurs first for interactions among the 
biggest bodies but later also between big and small bodies. The transition to RG is initiated and 



5.2.1. Delayed onset of runaway growth 



-45 - 



growth accelerates. 

For the 6 AU run the behavior can be regarded as a 'delayed' onset of runaway growth, i.e., 
RG takes off at Ri ~ 100 km but its general characteristics are not much different than the non-Ts 
run. However, the behavior of the 35VeGdFrTs model is more interesting. Here, RG also takes 
off at i?i ~ 100 km but then displays very extreme properties. Growth very rapidly produces a 
~2 X 10^ km oligarch (see Fig. Illb). The clue to the understanding of this behavior lies in the 
sizable number of fragments that are produced and in the different way these are treated in the 6 
and 35 AU runs. In the 6 (and 1) AU runs the fragments are assumed to move with the gas at 
a fixed approach velocity Va = riv^. However, at 35 AU such a restriction was not applied, i.e., 
fragments decouple from the gas and cool themselves efficiently through inelastic collisions. Thus, 
we have the peculiar situation that the random velocity of the big bodies [vm) is larger than that 
of the smaller bodies (here: fragments), Vf. 

At the point where the massive particle fulfills the condition vm < Vesc,M dynamical friction 
with the low-mass bodies (and fragments) starts to further decrease its random motions vm- Due to 
the enhanced GF, these bodies quickly accrete the fragments. The situation is exacerbated because 
the relative velocity between the fragments and the runaway body is set by vm (i-e., the random 
velocity of the runaway body), and not Vf (the random velocity of the fragments), since vpi > Vf. 
Since vm is decreasing, the growth displays more extreme characteristics with RG-index k > 4/3 
(see Sect. 113]). 

Note again that for the 35 AU Ts run we have treated an academic case since turbulent stirring 
by density fluctuations in the gas disk cannot operate on timescales longer than the disk dissipation 
timescale. Another word of caution concerns the (artificial) sharpness of the transition between 
the superescape and dispersion-dominated regime; in our model dynamical friction is suddenly 
'switched on' at v = Vesc, whereas in reality the transition occurs smoothly. Notwithstanding 
these concerns, this simulation can be regarded as representative for runs that are characterized by 
initially superescape velocities, but which are nonetheless accretionary. Accretion timescales are 
generally long; however, at the point where vm ^ Vesc,M RG sets in and a very rapid evolution 
follows, due to sweep-up of a population of dynamically cold fragments that has been produced 
during the preceding superescape phase. 

5.3. 35 AU gas-free models 

We analyze the 35VeFr run, which includes fragmentation but no turbulent stirring, in some 
more detail. Figure [T3k provides the mass spectrum at several times. During the runaway phase the 
high mass tail of the column density distribution is characterized by a power-law slope with index 
p < —2. This slope fiattens during the RG phase and, like in Fig. [8l approaches a value p —2.5. 
However, the evolution top = —2.5 does not fully complete at the higher masses; compared to Fig. [8] 
it seems to break at an earlier stage. We attribute this difference to the presence of fragments that 
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Fig. 13. — (A) The surface density distribution at eight times during the runaway growth and 
ohgarchy stages of the 35VeFr run. The fragment distribution is not shown here. (B) The velocity 
distribution for the same times as in (A). The dashed hne indicates the slope for dynamical equi- 
librium, while the thin solid lines gives Vh{R) = RQ/a and fosc(^)- The times are indicated in the 
legend. For the t = 7.7 x 10^ yr curve two auxiliary lines are plotted, illustrating (1) the position 
of Ri{t) and (2) Vh{Ri) at this point. The most massive particle accretes bodies to the right of the 
dot in the shear-dominated regime. 



start to dominate the accretion behavior (see below). 

In Fig. 113b the velocity distribution is plotted. It can be seen that initially the distribution is 
in dynamical equilibrium, for which v oc oc (dashed line). However, towards the end of 

the simulation the low-m bodies no longer obey this relation; the velocity distribution flattens out 
at low R. The reason for this behavior is that these bodies are stirred faster by the bigger bodies 
than they can equilibrate by dynamical friction with other small bodies. The random velocities 
among the small bodies lie close to the escape speed, which implies that GF and subsequently the 
interaction rates are weak. Dynamical friction among these bodies is therefore suppressed. On the 
other hand, the stirring by the runaway bodies/oligarchs does not discriminate between the mass 
of the small bodies: the big bodies regard all bodies at lower m as (massless) test particles. 

The lower, thin solid line gives the Hill velocity, Vh{R)- Bodies of size R that lie below this 
curve can accrete other, less massive bodies in the s.d. -regime. These are mostly particles of similar 
mass, but for the most massive body the s.d. -regime applies for almost one order of magnitude in 
size. Similarly, fragments (not shown in Fig. 113b) are mostly accreted in the s.d. -regime, because 
collisions among fragments keep their random velocity low. 

Figure O presents the spatial distribution of all particle groups the simulation contains at 
three different times. See the caption of Fig. [9] for the coding of the symbols and colors. The 
decreasing v/vh during the runaway growth stage can clearly be seen from the bluer colors in 
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Fig. 14. — Scatter plot of bodies' mass and position at tliree times during tlie 35VeFr simulation. 
See the caption of Fig. [Uj for the description of the symbol- and color-coding. 

Fig. 114b compared to Fig. [T^ . However, the random velocity at a given R always increases with 
time (Fig. 113b) and collisions among low-mass bodies result in fragmentation since v exceeds the 
escape velocity of these bodies. Runaway growth is clearly fast {vh grows faster than v) but it 
is also clear that strong particle separation - a key signature of RG - does not take place due to 
the fact that there are so many competing bodies. In Fig. 114b . and especially in Fig. [T^ . this 
separation is more obvious, however. But by now the system is in oligarchy: a few bodies have 
separated from the main distribution and these are dynamically heating the remainder. The colors 
turn red again. 

5.3.1. Which bodies contribute to the growth? 

It is instructive to see what kind of collisions dominate the growth of the most massive particles 
during the RG and oligarchy stages. This information is presented in Fig. [15] at three intervals, 
during which the size of the most massive particle, doubles. For example, the lower (light grey) 
symbols give the contribution by mass of the particles that accrete with the largest body over its 
growth from Ri = 60 km to Ri = 120 km. The detached square gives the contribution from the 
mm-size fragments. What can be seen from Fig. [TJ] is that during the RG stages {grey curves) 
bodies at each R < Ri contribute approximately equally to the growth of the largest body. The 
assumption that only the bodies that dominate the mass (i.e., those at i? ~ Rq) contribute, would 
be wrong. Instead, collisions that take place in the shear-dominated regime, i.e., those in the tail 
of the distribution, contribute a sizable fraction to the growth of Ri, despite their low abundance. 

In the final stages {black curves) this trend reverses. Here, bodies at ~i?o and contribute 
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Fig. 15. — The contribution by mass of bodies accreted by the maximum particle. The dotted lines 
indicate the radius-mass relation for a factor of 10 increase in the mass of the most massive particle. 
The corresponding starting and end points of Ri{t) are indicated by circles. For each case, the 
distribution of collision partners during the growth of Ri{t) are shown (crosses). The contribution 
of fragments is shown by squares. The data are averaged over 5 simulations. 



most with the contribution from intermediate-size bodies being suppressed. This signifies that the 
transition to the (classical) two component oligarchy state has taken place. In fact, the contribution 
from fragments (indicated by squares in Fig. \TE\i starts to dominate. But the key insight is that 
during the RG-phase the two-component approximation is invalid: all mass ranges contribute to the 
growth. As explained in Sect. 13.31 and Sect. 14.21 scattering among the high-mass bodies should not 
seriously affect the validity of these conclusions as long as the density of oligarchs is high enough. 



Makino et al.l (j 19981 ) recognized that the large focusing factors among particles in the high 
mass tail of the size distribution compensates for their low numbers. Under the assumptions that 
the relevant qu antities, i.e., th e veloc ity spectrum and the collision radii, can be given as power-laws 
of their masses, iMakino et al.l (| 19981 ) solved for the steady-state value of the mass-distribution and 
showed that it was consistent with a power-law index of p = —8/3. This is consistent with what 
we find, although we would like to emphasize the dynamic nature of the process. During RG, a 
power-law at the high ma ss tail of the distribution emerges until a point is reached at which it 
breaks fcf. IWetherilill99ol ). 



In the Makino et al. (jl998l ) model collisions among the largest bodies contribute, most to 
the growth of the most massive body. This is in contrast to our results, where we find that all 
bodies contribute roughly equa ll y. Th e reason for this discrepancy can be found in the power-law 
assumptions that Makino et al. (jl998l ) employ, i.e., that thermal equilibrium holds at all sizes and, 
more critically, the neglect of the shear-dominated regime. Under these (idealized) conditions, 
GF-factors become infinite, which in reality is not possible. 
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Fig. 16. — Outputs of various simulations without gas drag but with fragmentation, performed at 
35 AU. 

5.4. Miscellaneous models 

To test the robustness of our simulation results against the initial conditions (random velocity, 
planetesimal size, surface density) we performed additional runs at 35 AU in which these quantities 
(S, Rq and vq) are varied. The results are presented in Fig. [161 The y-axis in Fig. [16] is again 
normalized to trun (Eq. [23]) which is a fixed but different number for each run. 

Perhaps the most striking feature is that the normalized timescales are all of order unity for 
all models, indicating that trun (or Tj-g) is the appropriate timescale to characterize the runaway 
growth phase. Furthermore, the shapes of the curves are generally similar: all runs show a peak in 
-M^i/n^* (Fig. [T6k ) and a minimum of Vx/vh (Fig. [T6b). The associated steepening of the timescale 
curves is not always so clear in Fig. [T6b : growth becomes fragment-dominated in many simulations 
and remains fast. 
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The radius where the minimum of Vx/vh takes place varies. Compared to the 35VeFr curve it 
seems to occur at a lower size (~500 km) for the low Rq and low S runs but at a larger size (really 
towards the end of the simulation) for the high Rq and high E runs. For example, comparing the 
low and high S runs in Fig. [T6b (dotted curves) one sees that they diverge at around Ri ~ 300 km. 
The high-S run stays longer in the RG phase. 

Regarding the variation in initial velocity vq , note that the Iow-uq run (dashed grey line) quickly 
catches up with the standard run (the black solid line in Fig. 116b). The reason is that planetesimal- 
planetesimal stirring is very effective for v < Vesc- This justifies our choice to determine the width 
of the zones from the escape velocity of the bodies (Eq. [5l Sect. 12. 3p . 

The high velocity run (black dashed curve), on the other hand, shows very different behavior. 
Initially it is (expectedly) much slower than any other model but then at i?i ~ 20 km there is 
a sudden transition to a much faster growth mode. The behavior of this run falls in the same 
category as the 35 AU turbulent-stirring one, explained in Sect. 15.2.11 The superescape velocities 
produce copious amounts of small fragments and once the big body is able to use its GF to sweep 
them up, (extreme) runaway growth follows. 



6. Discussion 



6.1. Transition from the runaway growth to the oUgarchy phase 



In Sect, m we have identified the minimum of Vx/vh, the point where gravitational focusing 
(GF) peaks, as the transition between the runaway growth (RG) and oligarchy accretion phases. 
The corresponding transition size, Rtr, is also found to be close to the point where Mi/m^, peaks. 
We have tabulated i?tr for a variety of simulations, corresponding to varying physical conditions at 
several disk radii (Table [5]) . 



We now compare i?tr with the prediction by llda and Makind (119931 ). Eq. [H Using Sm = 
M/(27raAast), ~ S the density in solids, with Acgt = ARh the width of the heating region 
where A k, 5 reflects the spacin g among oligarchs, one obtains for the RG/oligarchy transition of 
Eq. n (cf. iThommes et al.ll2003l ) 



R 



Tg/oli 



^AaT.Rl 



94 km 



g cm 



_2_ 

" 15 



10 g cm 



-2 



f— V 

VAu; 



Rn 



10 km 



(24) 



Comparing the theoretical prediction Eq. [25] with the 'measured' transition points from our 
'indicators' (i?* and i^tr, see Table [5|) we see a clear discrepancy; typically Rtr is larger than Eq. [H] 
by several factors in radius but the discrepancy increases for the outer disk models. For example, 
whereas according to Eq. [23] RG will stall at the ~100-200 km size, we see from Table [5] that i?tr 
is rather ~300 km for the 1 AU models, ~600 km for the 6 AU models, and '^lO^ km for 35 AU 
models. Clearly, the domain of pure RG extends a little further than Eq. [2^ predicts, especially in 
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Fig. 17. — The runaway growth timescale Tj-g normahzed to trun = RoPsf^^ for various runs hsted 
in Table m The standard runs including gas drag and fragmentation (first 9 rows of Tabled]) are 
indicated by circles. Runs characterized by large random velocities are indicated by diamonds 
and other 35 AU runs are indicated by crosses. The majority of the runs are well represented by 
Kj-g ~ 0.03-0.1. But the high velocity models deviate significantly from this trend. 



the outer disk. 

The underlying reason for this discrepancy is the fact that RG is irreconcilable with the two 
component picture Eq. [24] relies on. For example, we have seen in Sect. I5.3.1l that the biggest body 
accretes from all mass bins, not just from the ones that dominate the mass of the distribution. 
The same holds for the stirring of small bodies. For a p ~ —2.5 mass spectrum the biggest bodies 
dominate the stirring, but it is not one big body that dominates. In other words, Eq. which 
presumes that the two component approximation is valid, cannot be applied to the RG stage. 

However, for oligarchy a two component approximation becomes valid. Indeed, it may be 
defined as such. Therefore, the start of oligarchy can be defined at the point where the relevant 
timescales in the two component approximation match the RG timescale Tj-g. Initially, Tj-g is the 
dominant (shortest) timescale and the two component picture is invalid. However, due to the 
increasing GF-factors there will be a 'tipping point' after which a two-c omponent picture does 



become applicable. We have recently addressed this issue quantitatively (jOrmel et al.ll2010l ) and 
found an analytic estimate for the transition size, -Rtr^ 



0.1 y Vigcm-3y Viokmy vau/ Viogcm-2 

Using the values for S and Rq listed in Table H] Eq. [2^ gives transition radii of Rtj- ~300 km, 600 
km, and ~10^ km for the runs at 1, 6, and 35 AU, respectively. This corresponds reasonably well 
with the obtained transition radii in Table [5] Note the rather weak dependence of Rtr on Kj-g. 
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In Eq. [25] K^a- is related to Ti-g as in K^g = T^g/t^un where t^un has been defined in Eq. [23l In 
a dimensionless form, K^g characterizes the outcome of our RG simulations. In Fig. [T7]we plot its 
value for many of the simulation runs listed in Table [5j Overall, we find that K^g lies in the range of 
0.03-0.1 for simulations that include the relevant physics. Since K-^g is reasonably well constrained, 
Eq. [25] is a robust prediction. However, simulations that are characterized with very large (initial) 
random velocities or include turbulent stirring do not obey this trend. But these simulations are 
(initially) not in RG. 



6.2. 



The importance of scattering and migration of bodies (neglected effects) 



Scattering and migration have been ignored in this study. In Sects. [3^3] and [^121 we have already 
discussed the effect of gravitational scattering of bodies. In our simulations we treat only a small 
patch of the disk and we do not expect scattering to significantly change the outcome. However, 
the assumption that a local representation of the disk is applicable itself may be questionable. 
Scattering of bodies or long-range interactions from the inner disk (where evolution timescales are 
shorter) may affect the evolution of bodies in the outer disk. In future studies, we may include the 
change in semi-major ax is due to the scattering of bo dies in our model, following, for example, the 
prescriptions outlined in [Tanaka and Idal (Il996l . 119971 ) . 



Similarly, we have neglected migration of big bodies and subse quent shepherding o f small 
bodies. We can, however, assess the timescale due to type-I migration (jTanaka et al.ll2002l ): 
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for a solar mass star, and a protoplanet size i?pp of internal density = 3 g cm-3. Since R 



■-pp 



2 X 10 km is the maximum size we reach in our runs, the effects of type-I migration are rather minor 
as Tjnigr is much larger than any relevant timescale we consider. However, for the subsequent core- 
accretion p hase, where -Rpp gr ows to Earth-size proportions or larger, type-I migration becomes 
important ([Tanaka et al.[ [20021 ). 



Gas drag-induced migration is also not included in our manuscript. The timescale for migra- 
tion of km-size, planetesimal bodies, may still be reassuringly long not to affect any of our key 
conclusions, but for small fragments it is another story. Th ese could be rernoved f rom the region 
where they are produced in as little as a few hundred years (jWeidenschillingl [l977a[ ) . We will now 
take a more closer look to the behavior of the fragments. 



6.3. The importance of fragmentation 



In simulations where a lot of fragments are produced, like in the runs that include turbu- 
lent stirring or a large initial velocity dispersion, the shear-dominated (s.d.) regime can become 
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important. In Appendix IB.1.31 we obtained dM/dt and , using Eg. [201 w e find a timescale (cf. 
Goldreicli etaPbooi iRafikovl booi IWeidensdiillinal bood : Ichamberslboog ) 
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(27) 



Equation [27l assumes that the fragments constitute a cold, very thin disk {v < Vh and Vz < a^^^). 
The striking feature of Eq. [27] is that, except for S, it does not depend on disk radius o. In the 
outer disk, s.d. -growth by sweepup of fragments is the only fast growth mode available. Despite 
the scaling with i?pp, which shows that growth is not in the runaway mode, it is here much faster 
than its runaway counterpart, i.e., TI^^^'"^- <^ T^g. 

However, there are many caveats regarding our treatment of fragments. First, there is the 
uncertainty in the collision model, reflected in the parameters Ofr and e (see Sect. ISTT]) : changing 
these parameters may increase or decrease the significance of fragmentation (and its re-accretion). 
Another concern is the presence of gas drag. Equation 1271 assumes that the fragments, like the 
protoplanet, move at the Keplerian orbital velocity. We have already mentioned the fast radial 
orbital decay these particles experience. Furthermore, if the particles are small enough they move 
with the gas at a relative motion of rjVk ~ 30 m s~^ with respect to the planetesimal bodies. 
Interactions between massive bodies and fragments then probably take place in the high-velocity 
regime ('dispersion-dominated' is perhaps unsuited here since the velocity difference is systematic) 
and Eq. [27] is inapplicable. In our simulations, we have tried to pre-empt this effect by fixing the 
random velocities of the 1 and 6 AU runs at v = ijVk- We are currently investigating in more 
detail the effects of gas drag on the gravitational cross sections for small particles (Ormel & Klahr, 
submitted). 

It has previously been reported that the settling /dampi ng process of gas drag allow s for a very 
efficient growth of fragments in the s.d. -regime via Eq. [27] (jKenvon and Bromlevl l2009l ) . However, 
these fragments should occupy a very special niche: they cannot be too small as they would 
otherwise couple too strongly to the gas to allow efficient accretion, but cannot be too big either 
since then the assumption of a thin disk and efficient cooling becomes problematic. Thus, it is 
unclear whether Eq. [27] can materialize in gas-rich environments. On the other hand, it is likely 
that the coUisional cross section Rm] is enhanced above the 2-body limi t due to the dissipative 
natu re of the encounter (llnaba and Ikoma] l2003l : iMuto and Inutsukal [2009j; iTanigawa and Ohtsuki 
2010l ). This is an active area of research. 



7. Summary and Conclusions 

We have developed a new statistical code to study the evolution of the planetesimal size 
distribution. The key novel element is that the code treats interactions between particles, rather 
than mass bins. We have tested the code against existing literature A^-body and statistical studies. 
Starting from an initially monodisperse distribution of ~km-size planetesimals, we have performed 
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a large parameter study with the aim to identify the transition between the initial runaway growth 
phase and the later oligarchy phase and to assess the sensitivity of key physical processes like 
gas drag, fragmentation, and turbulent stirring on the planetesimal accretion process. Our main 
conclusions are the following: 

1. Interactions that take place in the dispersion dominated regime, where random velocities v 
lie between the Hill velocity Vh and the escape velocity Vesc of the largest body, fulfill the 
conditions for runaway growth, Eq. [T9j However, for spatially resolved systems, one must 
distinguish between the local and the global variant of the runaway growth definition. In 
oligarchy, the system is only locally in runaway growth, i.e., Eq. [19] no longer holds between 
bodies that are spatially isolated. 

2. In simulations that start out at the very typical condition where the initial random velocity 
vq of the bodies is less then their escape velocity, fesc,0) runaway growth ensues. We find 
that during the runaway growth phase all masses contribute - approximately equally - to the 
growth of the largest body, which occurs exponentially at a characteristic timescale Trg. 

3. The runaway growth timescale Tj-g is empirically determined and expressed in terms of the 
initial parameters as T^g = Kj-gPsRo/T,Q. We find a value of Kj-g of 0.03-0.1 is typical. 
Simulations that are characterized by high initial velocities {vq ^ Vesc,o) do not obey this 
trend. 

4. We find that at a transition size -Rtr the runaway (exponential) growth is over and becomes 
much slower (a considerable distribution of fragments may mitigate this effect, though). At 
this transition point the gravitational focusing factors peak. 

5. During the runaway growth phase the mass distribution at the high mass end gradually 
changes into a power-law Ns(m) oc with p approaching ~— 2.5 near the end of the RG 
phase. During the oligarchic phase the distribution breaks and becomes characterized by two 
components: the oligarchs (at high m) and leftover planetesimals. 

6. Due to the gravitational stirring, interactions among low-mass bodies reach the fragmentation 
regime v > Vesc- In this study we have treated such collisions as erosive, producing copious 
amounts of ~mm-size fragments. The ability of fragmentation to influence the growth of 
the biggest bodies (through re-accretion of fragments) depends somewhat on the adopted 
collision parameters. Fragments do not dominate the evolution during the runaway growth 
phase, except in models that include (signiflcant amounts of) external stirring. However, 
during the oligarchic phase collisions become violent enough (and timescales long enough) for 
fragments to become important. 

7. Simulations that start out at v > v^sc but where velocities are still sufficiently low to be 
accretionary are characterized by long accretion timescales and copious amounts of fragment 
production. At the point where the escape velocity of the largest bodies starts to exceed 
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Table 6: Summary of runaway growth scenarios 
Case Description 

(A) Classical regime 

Prerequisites: Velocity dispersion below escape velocity of biggest body, {v < v^sc) 

Key characteristics: Column density evolves into a power-law distribution, Ns{m) oc 
with p ~ —2.5 

Collisions between all masses contribute to the growth of the biggest 
body 

Runaway growth fast with timescale given by Trg = i^rg^run (see Eq. [23|) 
with Krg ~ 0.03-0.1 
Outcome: Transition to oligarchy at radius given by Eq. [25] 

(B) Fragmentation-dominated 

Prerequisites: Production of a sizable amount of dynamically cold fragments with vj < 

Key characteristics: Very fast growth possible (shear-dominated regime) at timescale given 
by Eq. [27] 

Outcome: (most likely) 2 component oligarchy of protoplanets and fragments 

(C) Superescape regime 

Prerequisites: Velocity dispersion above escape velocity of the biggest body {v > Vesc), 

but net accretion 

Key characteristics: Continuous size distribution, declining exponentially at high-m; no grav- 
itational focusing, slow growth, fragmentation 

Outcome: Transition to runaway growth (scenario A) at point where v < Vesc', 

possibly strong fragment-dominated growth (scenario B) 



V the growth mode turns to runaway. It becomes especially fast if the fragments are kept 
dynamically cold due to mutual (elastic) collisions. 

8. Sweepup of such a dynamically cold population of fragments takes place in the shear-dominated 
regime, which can lead to very rapid growth rates (see Sect. 16. 3p . However, we have ques- 
tioned the viability of this mechanism in gas-rich systems since particles are tied to the gas 
and suffer radial orbital decay. 

From these general findings, we construct three scenarios through which planetesimal growth 
could have proceeded, see Table [6] In the first scenario, the classical regime, the system starts out 
in the dispersion-dominated regime, Vh < v < Vesc, and runaway growth ensues according to points 
2-5 listed above. At Ri{t) = i?tr (Eq. [25]) the runaway growth phase is over and is superseded by 
oligarchic growth. This is the point where a full 2-component approximation (oligarchs and smaller 
bodies) becomes first applicable. Most studies that deal with planetesimal accretion have focused 
on this scenario. 
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Planetary accretion may have deviated from the above-sketched contours of the 'classical 
regime'. For example, if a populous reservoir of dynamically cold particles (fragments) is present, 
the accretion takes place in the shear-dominated regime {vx < Vh, scenario B). In the shear- 
dominated regime gravitational focusing factors are not determined by the random velocity of the 
particles, i.e., the self-regulated effect which is the hallmark of oligarchic growth is absent. We ob- 
tain such a situati on if the initial condit ions features an already mature body embedded in a 'sea' 



of small particles (IWeidenschillingi l2005l ) or when collisions among the (rubble-pile) planetesimals 
lead to fragmentation (i.e., v > Uesc)- Provided the fragments are able to cool themselves (through 
mutual highly dissipative inelastic collisions) they could quickly dominate the contribution to the 
mass gain of the biggest body. The shear-dominated regime marks a very fast growth mode, es- 
pecially for the outer solar system and for gas-free environments. However, fragments should be 
kept dynamically cold and we have raised a note of caution on the viability of this mechanism in 
gas-rich environments. 

The final scenario (C) is characterized by random velocities of the planetesimal bodies that 
(initially) exceed the escape velocity {vx > fesc) but are nonetheless accretionary. External stirring, 
e.g., as a result of density inhomogeneities in the gas disk, may provide these conditions (provided 
it is not too strong to shut-off accretion altogether). Since gravitational focusing is negligible when 
V > fesc) accretion timescales are very long. As the system is not in runaway growth, the size 
distribution remains continuous with an exponentially-declining tail at the large masses, rather 
than a power-law. However, at the point where the largest body fulfills Vcsc,Ai > vm, growth will 
enter the runaway regime (dispersion-dominated regime). If a lot of fragments has been produced 
and mutual elastic collisions are able to keep the random velocities of these particles low, the 
runaway effect is especially pronounced since the gravitational focusing factors are determined by 
the random velocity of the biggest body {vm)- Eventually, growth could enter the shear-dominated 
regime (scenario B). 
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Fig. 18. — Sketch illustrating how the zoom factors are determined, {left) The particle distribution 
of a certain zone at time t is binned in terms of log mass. The group size Ng and the zoom factor 
z of the swarms are determined from the number of particles in the bin Np^y^^f, and the particle 
resolution A'res- In this example, A^res = 80 and 8 RBs are assigned per mass bin. For the z = 
bins the high-mass bodies are individually resolved, (right) Illustration of the state in the indicated 
bin, for which the group size is 4 (z = 2, dashed rectangle). One has that qa + 93 + gc = -^physj 
with gi the true number of particles for a species group i. Particle swarms A and B are resolved, 
C is not. 
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A. Determination of N* 



In Sect. [2^ we introduced Ng, the number of physical particles that are associated to a rep- 
resentative body (RB). The group nu mber, A''^, is deterinined by a function that depends on the 
distribution (the distribution method: lOrmel and Spaansll2008l ). As the distribution changes with 
time the amount of grouping likewise adjust. We denote this function by N*. Here, we sketch how 
it is obtained. 



This preprint was prepared with the AAS lATf^X macros v5.2. 
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We illustrate the distribution method using Fig. [18] as an example (the numbers given here 
are entirely arbitrary). The distribution method assigns the RBs equally over logarithmic mass. 
Therefore, we first bin the particles by mass. This binning is just an auxiliary feature and not 
fundamental to the program, i.e., it is only used here for the determination of A'^. The zone that 
Fig. [T8] corresponds to is covered by in total 80 RBs (A'^rcs = ^rh/^zo = 80). These are equally 
assigned over the 10 logarithmic mass bin, i.e., 8 RBs per bin. Therefore, the size of the groups 
equals Ng = A'^phys.fc/S where A^phys.fe is the number of physical particles in mass bin b. Numbers 
are rounded where needed; i n fact, Ng is taken to be a multiple of 2, i.e., Ng = 2^ with z integer 



(see lOrmel and Spaansll2008l for rationale). Thus, the result is that the numerous low-mass bodies 
are highly grouped (large Ng), whereas the high-mass bodies are individually resolved {Ng = 1 or 
z = 0). This property allows RG to be accurately modeled. 

In the right panel of Fig. [18] we further illustrate the relation between the RBs and physical 
particles. We focus on bin #7, which has Np^j^j = 32 physical particles that are represented by 8 
RBs of Ng = 4. In principle, each physical body is different from another, i.e., it is characterized by 
unique physical properties (mass, velocities, etc.). However, the program only deals with encounters 
between RBs - i.e., particle groups that share the same properties - which preserves an inherent 
level of 'graininess' to the distribution. We will refer to these ensembles of identical particles as 
'species'; e.g., there are qa physical particles of species A that have identical properties which are 
different from those of species B. 

Thus, we see that since the group size equals iVg = 4, 5 RBs are assigned to the particles of 
species A, 2.5 to those of species B, and 0.5 to those of species C. The total amount of RBs add up 
to 8, while the total amount of physical particles add up to Npi\ys,7 = 32. The fact that the total 
number of RBs assigned to the B-particles is not an integer does not pose a problem, as long as Ng 
particles can take part in the group collision. However, for the 'C-particles' the number of physical 
particles falls below the group size. These situations are strictly forbidden; under resolved groups 
are merged with other swarms which are closest in mass and velocity space (and also in radial 
position: merging occurs only for groups within the same zone) according to the criteria discussed 
in Sect. 12.21 For example, the C-group may be merged with the B-group, provided their properties 
(masses, positions, eccentricities, etc.) do not greatly differ. 



B. Calibration of the model 

The interaction timescale of a single body with a group of bodies j is tint = ("-jCint^^a)"^, 
where rij is the number density of j-bodies, fiint the interaction cross section, and Va the approach 
velocity. Inverting this expression, the interaction rate is defined as 
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where we have substituted dint = t^RxRz and Uj = Nsj/2h(,ff. In Eq. IBll Ngj is the column 
density of j-bodies, h^s the effective scaleheight, and Rx, Rz, respectively, the planar and vertical 
interaction radii. These latter quantities follow from R^^t md from the geometry of the encounter. 
For, example, Rz can never exceed the scaleheight hes, i.e., Rz = min(/iefr, i?int)- For each type of 
interaction, furthermore, the interaction radii depend on the velocity regime. The superscript '1' 
in A^^^ emphasizes that Eq. IBll is the collision rate for a single particle. 

Using Eq. IBll accretion rates, dM/dt, and stirring rates, dv'^/dt, can be constructed. For 
accretion we simply multiply Eq. IBll by rrij , while for stirring we multiply by the change in v"^ that 
the body suffers during the encounter, see Appendix IB. 21 i.e.. 



dv^ -KRxRzVa^j 2n 



int ) 



(B2) 



where 'int' refers to either dynamical friction or viscous stirring. Using these expressions and the 
procedure sketched in Sect. 12. 4^ we will verify that the geometrical model is, within factors of unity, 
consistent with more refined literature treatments regarding accretion and stirring. However, con- 
cerning the numerical application of the model we do regard these offsets as important. Retrieving 
these order unity factors by comparing with existing literature expressions is what we understand 
under 'calibration' and is the purpose of this section. 

A swarm of bodies of individual mass m is often characterized by a distribution in random 
velocities. In particular, the Rayleigh distribution is frequently adopted; i.e.. 



4v'v'^ 



exp 



(B3) 



(jGreenzweig and Lissaueijll992l ) is the probability that a body has random velocities given 
that V and Vz are the rms-values of the distribution. One would naturally expect that the mutual 
interactions within a population produces a distribution, a nd A^'-body experiments indeed indicate 



that this is the case ( Wetherill 



198C 



Ida and Makindll992l ). For this reason, accretion and stirring 



rates are often given as distribution-averaged quantities. When calibrating our model we will follow 
this convention - i.e., calibrate against the distribution average. 



B.l. Collision rates 



B.1.1. The superescape regime, w > Ve 



In this regime we assume (3 = i/e = Vz/v = 0.5 (jida et al.lll993l : iTanaka and Idalll996l ) and 
that the scaleheight of the disk hes exceeds Rcoi such that Rz = Rx = Rco\ = Rs- Furthermore, 
Va = V and Va/heg = v/{vz/^) = ^/I3 = 2Q and the accretion rate Eq. IBll becomes 



(B4) 
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where we have augmented Eq. IBll by Ai, the cahbration factor. We c an compare Eg. IB4I with 
the one-bod y accretion rate, averaged over the distribution, obtained bv iGreenzweig and Lissauer 
(|l990l . I1992I ) 



(A 



R 



N. 



1 ' 2tt 

where J- {(3) is an integral expression of f3. A numerical evaluation gives J-'(0.5) 
we find Ai = 0.90. 



(B5) 

16.1, from which 



B.1.2. The dispersion- dominated regime, 2.5vh < w < Vesc 
In the dispersion-dominated regime, i?coi = RsVesc/w. For our calibration we again assume 



Rcoi = Rx = Rz = Rs < hes and Va/hes = 2r2. Then, 



dt 
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with A i the calibration factor for this regime. The equivalent expression bv lGreenzweig and Lissauer 
\imi ). their Eq. 17, reads 



^dNacc^ 



R^M 



w 



dt ' 2ir 

A further numerical evaluation gives G{0.5) = 43.0. Ignoring the J-' term, then gives A2 = 1.5 



(B7) 



B.1.3. The Hill regime, w < Vh 



In the Hill regime, the random relative velocity w of the particles falls below v^', the approach 
velocity is now determined by the Keplerian shear, Also, for the calibration we assume 

that the vertical accretion radius is limited by the disk scaleheight, Rz = /leflf- Not all particles 
approaching at impact parameters |6| < = 2.5Rh manage to enter the Hill sphere, however. 
Particles approaching at small impact parameters move on horseshoe orbits and their trajectories 
will move away from the Hill sphere. Numerical studies have shown that only particles approaching 
at impact parame ters 1.7 Rh < b < 2.5Rh, i.e., only a fraction fin ^ 0.3 of Rrp], wi l l make it into the 



Hill sphere (e.g., iNishidal Il983l : IPetit and HenorJ Il986l : llda and Nakazawal 1 19891 : iGreenberg et al 



1991 



). The average approach velocity of these particles is ^flb = •i.2vh and we take this value as 
the approach velocity in the Hill regime (see Fig. [3D. Finally, only a fraction /hit of the particles 
that enters the Hill sphere will be accreted. Combining these expressions, we obtain 



A 



(1) 



/in/hit^: 



irRxRzVaN, 



2h, 



ef 



2,.77RhVhN,jY.^. 



(B8) 



We estimate /hit rather crudely by assuming that within the Hill sphere the particle motion is 
random at an average velocity of ~ 2.5vfi. This translates to a 2-body impact parameter of 
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Fig. 19. — The colhsion rate X^^^ /Nsj expressed in Hih units as func t ion of Hih inclination, 
ih = Vz/vh- Symbols show the numerical results of llda and Nakazawal (|l989l ). where we have 
only included values ih,eh < 1, i.e., the low- velocity regime. For each ih, the eccentricity- values 
were averaged with the error bars indicating the scatter (which is small). The gray line shows 
Eq. IB8I with ^3 = 1 in /hit (Eq. IB9p . Setting the calibration factor to = 2.9 produces an 
excellent fit. 



6coi = RsVescf^-^Vh ~ a^^'^Rh- The hit probability is then simply the ratio of this factor to the two 
relevant length scales, 

/hit = ^3X^xminfl,^V (B9) 



where A3 is the three-body calibration factor. In the last term we have taken care that the hit rate 
is limited to the scaleheight of the particle layer (in which case the interaction becomes truly 2D). 



We compare Eq. IB9I with the findings of llda and Nakazawal (119891 ). llda and Nakazawal (119891 ) 
have numerically integrated trajectories of particles in the restricted three-body problem as function 
of inclination and eccentricity and obtained the collision rate for a = R/Rh = 10"^. Their results 
are presented by the symbols in Fig. [19] as function of the reduced (Hill) inclination, ih = v^/vh- 
In Hill units lengths are normalized to R^ and times to 17^^ such that the accretion rate Eq. 



becomes Xij^h = ^■^'^ fYiit^sj and Eq. |B9] /hit = A'ia^/'^m.m.{a^^'^ /ih,l). Since in the Hill regime 
the velocities are determined by the Keple rian shear, the out c ome is insensitive to e^', therefore, 
we have in Fig. [19] for each ih averaged the llda and Nakazawal (|l989l ) results over the eccentricity- 
values that have eh < 1- Figure [19] shows that when = 2.9 we obtain excellent agreement with 
the numerical results. Besides, our simple model correctly predicts the transition to the 2D regime 
when hcs < bcoi {ih < a^^'^)- 

In the case of a very thin disk, Vz/vh < a^/^, the ac cretion rate then becomes dM /dt 



V h.T,i, in agreement with pr evious estimates (e.g., iGoldreich et al.l l2004l : iRafikovl l2004l : 
Chamberal2006l : IWeidenschillingI l2005l ) . This 2D result is the fastest possible accretion rate in our 



model, but it does not qualify for runaway gr owth since the mass exponent, k (Sect. HTS]) . is less 
than unity, k = 2/3 (cf. iKokubo and Idalll996l ). 



B.2. Stirring rates: velocity change at interaction radii 

Before a comparison of Eq. IB2I with the literature can be made, we first have to specify the 
changes in velocity, Au^ and Ai;|,j, for dynamical friction and viscous stirring. This is the topic 
of this section. The results are summarized in Table [71 In Appendix IB. 31 we then perform the 
calibration. 



B.2.1. Dynamical friction (elastic collisions) 

We model the velocity change that results from dynamical friction with a ID fully elastic 
collision. The velocity after such a collision is given by 

, vm{M - m) + 2mn;^ mm ^ 

= — w^) — ' ^ 

= (mT^) • 

In Eq. IBIOI a bar reflects the orientation of the collision; e.g., = ivm, dependent on whether 
the collision is head-on (negative values) or tail-on (positive). Equation IBIOI leads to a change in 
the kinetic energy of 

MMvIj) = — — r^(l^ - vm){Mvm + mv;^); (Blla) 

[M + m)^ 

2 \ AmM 

i^{mv^) = — — (fAf - Vm){MvM + mvm)- (Bllb) 

{M + mj^ 

These expressions add up to 0, reflecting conservation of energy, regardless of whether the bars 
indicate positive values (tail-on collisions) or negative values (head-on collisions). In order to get a 
mean change we simply average these expressions over the head-on/tail-on collisions; 

A(t;l^)df = - mvl), (B12a) 

^{vDdi = ^M^^Y ^^'^M - mvl). (B12b) 
These expressions are used for the velocity changes that result from dynamical friction. 
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B.2.2. Inelastic collisions (bouncing or accretion) 

We optionally make the assumption that physical collisions fully dissipate their collision energy. 
When implemented, accretion and bouncing are modeled as a perfectly inelastic collision. For a 
fully inelastic collision we have that the velocity after the collision is the same for both particles 

VM = Vm= M + ^ ^ (B13) 

although in the case of bouncing they will not stick. Repeating a similar approach as above we end 
up with the velocity changes as 

A(t'if)ic = p7^^(-2^^M + mivl - vli)), (B14a) 
which shows that the total energy change M/S.{v\^) + mA(w^) is always negative. 



B.2.3. Viscous stirring: distinction between 2D and 3D interactions 

In free space, the absolute relative velocity before and after a two-body encounter is equal in 
magnitude. In that case, only dynamical friction operates. However, in the protoplanetary disks, 
due to the influence of the sun, the 2-body energy is not conserved. This holds even for encounters 
in the dispersion-dominated (d.d.-) regime, where the Keplerian potential acts as a reservoir with 
which energy can be exchanged. Prom this perspective encounters are in fact always three-body 
interactions and solar gravitational energy can be converted into random motion, and vice- versa on 
thermodynamic grounds. The energy exchange with the solar body preserves the total energy of 
the system, but not the energy of the two orbiting bodies as a subsystem. The process of extracting 
energy from (or adding to) the potential is better known as viscous stirring. 

To understand v iscous sti rring physically, we fo llow a geometrical argument that originated 



from ISafronovl (jl969l ) , see also iGoldreich et al.l (j2004l ) . We consider the case of a small body with 



random velocity Vm experiencing a close encounter with a more massive body moving on a circular 
orbit. Due to the encounter the phase angle of the smaller body will shift by, on average, 90 
degrees but the magnitude of the local relative velocity w will stay the same. The point is that 
this latter quantity is not equal to Vm'-, at quadrature it is, but if the encounter takes place at 
perihelion or aphelion the relative velocity of the elliptical orbit with the local Keplerian velocity 
\s w = I'm/2. Then, if the encounter takes place at one of the latter locations and changes the 
orbit towards quadrature, this will have circularized the orbit of the smaller body. For example, if 
i^a/p denotes the relative velocity at aphelion or perihelion and Wq that at quadrature, the above 
reasoning reads Wf^jp = Vm /2 = w'q = v'^ where primes denote the velocity after the encounter. 
Therefore, v'^ = However, if the encounter takes place at quadrature, could increase by 
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a factor two if the orbit is re-oriented towards perihelion or aphehon {vm = Wq = w'^i^ = 

and v'^ = 2v). On average, these reorientations result in a gain in random energy: the (absolute) 

change is larger at quadrature than at aphelion or perihelion. 

Prom this discussion it is clear that viscous stirring results from the (random) re-orientation 
of the phase angle of the bodies motion. We have defined the viscous stirring radius R^s such that 
Avm ~ Vm for the small particle. For the heavy particle, the encounter will lead to a response that 
is approximately a factor ~m/M less due to its larger inertia. More formally, we will weigh the 
change in velocity Av by the masses of the collision partners, i.e., an amount M/{M + m)vm goes 
to the small particle and an amount m/[m -\- M)vm to the big particle. Next, we recognize that 
these impulses lead to random changes in the phase angle and that only the squares add up, i.e., 
we define 

M 



m 



^^-m)^.= [-^^-) ■ (B15b) 

as the change in random energy resulting from a viscous stirring encounter. 

In the shear-dominated (s.d.-) regime, scatterings that take place within the Hill sphere can be 
strong; velocities of the small particles are boosted to ~Wh- Consequently, instead of the previous 
expression which holds for the d.d.-regime, we have (At;^)vs = {[M/ {M +m)\2.f)Vh)'^ and (Af|^)vs = 
([m/(m -|- M)]2.5wh)^- Here, we normalize to 2.bvh since this is the characteristic velocity for 
interaction in the s.d. -regime. 

In the next section, we will calibrate the resulting stirring rates against the stirring rates 



obtained by lOhtsuki et al.l ()2002l ) and invoke an order of unity correction, /vs- However, rather 
than merely a constant, we will see that /vs is different for horizontal {v) and vertical {vz) velocities 
and also depends on /3 = Vz/v. Viscous stirring is highly dependent on the geometry of the 
collision. In the s.d. -regime {w/vh ^ 1) interactions are typically 2D since the radius at which the 
interaction takes place (-Rvs) is typically larger than h^s- This means that eccentricities are more 
strongly excited than inclinations. The latter's amplitudes are correspondingly reduced by a factor 
{hes/Rvs)- Of course in a fully 2D setting one would not stir the inclinations. 

However, in the d.d.-regime, interactions are modeled as 3D. Viscous stirring then equalizes 
the random velocity c omponents, meaning that is being dr iven to an equilibrium value of /3 = 0.5 



for a Keplerian disk (llda et al.l Il993l : iTanaka and Idal 1 19961 ) . This can be thought of as a kind of 



equipartition, although as noted before there is no two-body energy conservation for viscous stirring. 
It is even possible for viscous stirring to produce negative rates when the ratio /? = Vz/v <^ 1, for 
which /vs will also become negative. 

We will therefore distinguish between the 2D and 3D regimes when next discussing the mag- 
nitude (and sign) of the viscous stirring factors /vs and /vs-z- 
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B.3. Stirring rates: literature calibration 



We compare our viscous stirrin g expression that follows from Eg . IB2la. n d tlie discussio n above 



against previous lit erature studies (llda et alj Il993l : iTanaka and Ida 



1996 



Stewart and Ida 2000; 



Ohtsuki et al.ll2002l ). For dynamical friction we will not introduce a calibration factor as we will 
find that the present formulation models the dynamical friction stirring rates reasonably well (see 
Appendix IB.4j) . 



According to lOhtsuki et al.1 (j2002l ) , the stirring of the m-particle in the limit of m ^ M is 
given byj^ 



dt 



(B17) 



where (-P)vs represents the dimensionless viscous stirring factor for eccentricity, averaged over the 
distribution. Similarly, (Q)vs encapsulates the stirring of inclinations due to viscous stirring. These 
functions are in turn functions of the (reduced) inclination and eccentricity, i/j and Vh- In addition, 
(P)vs and (Q)vs depend on the velocity regime, i.e., the shear-dominated {vm ^ Vh) and dispersion- 
dominated {vm ^ Vh) regimes. The s.d. -regime is modeled as 2D; interaction take place at a small 
angle, 6 = h^s/Rvs 1- The d.d.-regime is assumed 3D (0 ~ 1). 

We introduce the symbol /vs(/3) as the calibration constant for viscous stirring. Since (-P)vs 
and {Q)vs are different in the 2D and 3D regime we discuss the regimes separately. 



B.3.1. 2D regime (shear- dominated), h^s < Rv 



We assume the disk is flat compared to R^s, such that Rz = h^s- Inserting Ri^t 
Rz = hes, Va = 3.2t;/i, and Au^ = /vs(2.5t;/i)^ into Eq. \B2\ we obtain a stirring rate of 



dvl 



dt 



Jin NsMjvs{^-^Vh) = 23j^sNsMRhVh- 



2.5Rh, 



(B18) 



where /in = 0.3 is the fraction of particles that en ter the Hill sphere, see A ppendix IB. 1.31 lOhtsuki et al 
([2002) argue on basis of numerical integrations (IStewart and Idall2000l ) that (Pvs) = 73 in the s.d.- 
regime. Taking this value and equating with Eq. IB171 we arrive at a value of /vs ~ 3.0. 



For the vertical excitations, we have argued that fvs-z should contain a factor 6^ = (hpff/R 



{vz/2.bvh)^ = 0.16/3^(u/u/j)^, reflecting the geometry of the encounter. Similarly, lOhtsuki et al 



^ Here, we have rewritten Ohtsuki's expression in our notation, see Eq. 6 of lOhtsuki et al.l ( 20021 ) . For the 
eccentricity stirring, under the assumptions outlined above, we have (in Ohtsuki's notation) 

ri(ei)vs 



dt 



= alQNsjhi2{P)vs, 



(B16) 



where we have used that (m —)mi <^ 'nij{= M) with e\ — Vm/vk the eccentricity, ao the semi-major axis, Nsj the 
column density, and hi2 = (Rh/a)'^- Multiplying both sides by (a57)^ then gives Eq. IB17I 
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(j2002l ) gives (Qvs) ~ /vhY ■ Comparing these expressions we find 

RvsJ 

since then 23/vs-z equals {Qvs)- 



/vs-z = l.l(^V, (B19) 



B.3.2. 3D regime ( dispersion- dominated) , Vm > 2.5vh 

When Vm ^ Vh it can be assumed that Rz = Rx = Rvs ^ h^g and /3 = 0.5. Inserting 
Rl = mRh{vh/vraY, Va/2h^s = and At;^ = /vs(/3)?;^> Eq. [B2]becomes 

^ = 36vrAs(/3)«(^)'. (B20) 



Now lOhtsuki et al.l (|2002l ) gives for the high velocity regime 



vs 



dt TTChih 



logfl + A 



2n 



1 + A2 



(B21) 



where I^^iP) is an integral expression that depends only on /3 and A is the Coulomb factor. A similar 
relation exists for (Q)vs but then with The functions and I^s are plotted in Fig. [20l In 

evaluating the /vs (/?) terms we have approximated the formal definition (which involves an integral 
expression that cannot be analytically solved) by an exponential fit, /vs ~ oq + ai exp[— 02/?], see 
Fig. [200 

We approximate the term in the square brackets in Eq. IB21I as ~2 log A, which is appropriate 
if A 1. Inserting Eq. IB21l into Eq. IB 171 with e^ = Vm/vh and ih = Peh/vh = Vm./2vh then gives 



dvl 2m^,vlRlNsM 



dt TTV"^ 



{\og h)VL-^ (B22) 



and a similar term (but then with /^) for the vertical velocities. Comparing Eq. IB22l with Eq. IB20I 

gives 

( /vs \^^( lUP) \ 
\ /vs-z I I^siP) I ' 

where we have ignored the Coulomb term (see Appendix IB. 4p . 



(B23) 



Again, we see that our calibration factors are of order of unity. However, for low /3-values 
/^ becomes negative: eccentricities are strongly damped and inclinations strongly excited. The 
net effect of the negative eccentricity stirring is then to (rapidly) increase /3, until an equilibrium 



ChambersI ( 20061 ) applies a similar fit to the Pvs and Qvs expressions. 
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Table 7. Summary of the velocity changes that particles experience upon an interaction 



Case*^ Velocity change At;^ or Av 



M 



Collisions (Bouncing) Dynamical friction Viscous stirring*^ 



General Equation IB14I 

Vm--im<^ M) {vIj - vl,) 

Vm:{m = M) (-3f^ + z;|^)/4 

VM--{m<^ M) -2{m/M)vlj + {m/Mf 



Equation [BT2] Equation |Bl5] 

Avj^ - 4{m/M)vl vl 

Vm - vl, f^/4 

-4(m/M)f|f +4(m/M)2t!^ {m/Mfvl, 



'^For example, the special '■ [rn <C My means that the Au^ changes are given under 

the condition that m <^ M. 

^For viscous stirring, if the interaction takes place in the shear-dominated regime, t;^ has to be 
replaced by (2.5f/i)^. 



at 




Fig. 20.— The functions /^(/3) and as defined bv lOhtsuki et al.1 ([2002). Crosses denote 

numerical evaluations at discrete values of /?. The curves present a fit to IvJ^ of the form /(/?) = 
oo + CLi exp[— 02/?]. The fit parameters are indicated. 
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value /3* is reached, for which /3* = \J Ws/ I^s ~ 0.55. This behavior is quite the reverse from the 
s.d. -regime, where any inchnation-stirring is strongly suppressed with respect to the eccentricity. 
The key physical reason is that in the s.d. -regime the interactions take place at a radius much 
larger than the disk height, iZvs > ^effj which enforces the 2D geometry. In the d.d. -regime, on the 
other hand, i?vs < ^eff and interactions are 3D. Of course, our simplified treatment - to identify 
the s.d.-regime with the 2D case and the d.d.-regime with the 3D - does not fully do justice to 
the full complexity at the transition Vm ~ Vh- When a system moves from the s.d.-regime into the 
d.d.-regime, interactions are initially still 2D. However, the outcome will qualitatively be the same: 
i?vs dec reases, interactions b ecorn e increasingly 3D, and any stirring only enhances this trend. 
Indeed, iRafikov and SlepianI POIG ) have recently studied this specific setting and found that the 
2D d.d. -transition regime is extremely short-lived. 



B.4. 



Final expressions concerni ng st irring and comparison vi^ith lOhtsuki et al 

(j2002l ) stirring curves 



We summarize the expressions for stirring rates that result from the geometrical model. These 
can be concisely written as 




ttR^RzN, 



/vs(/3) 
/vs-z(/3) 



(Az;2)i^tlogA. 



(B24) 



In the s.d.-regime /vs = 3.0 and /vs-z is given by Eq. IB19[ In the d.d.-regime /vs(/3) is given by 
Eq. IB23I Furthermore, the d.d.-regime is characterized by a Coulomb term, log A. The Coulomb 
term takes account of the encounters that occur at Rint < h < hes, which collectively can give 
a contribution - which is of order unity, but important for the numerical implementation of our 
work. We estimate the magnitude of the contribution to be determined by the ratio of R^^t to the 
scaleheight. 



log A = log exp[l] + 



Rh 



(B25) 



(the inclusion of the exp[l] term ensures that log A > 1 for all values of /leg /-Rint) 



These expressions can be compared with other works, e.g., lOhtsuki et al.l (j2002l ). Figure [21] 
presents the stirring rates that follow form Eq. IB24l for f3 = 0.5. Rates are given in Hill units, where 
lengths are normalized to Rh and times to il"^ and A'^^j = 1, see the discussion in A ppendix IB. 1.31 



One t hen arrives at the dimensionless expressions denoted by P^s, -Pdf > stc, for which lOhtsuki et al, 
(j2002l ) gives analytical fits. These are plotted in Fig. [211 



Comparing the curves, one observes a satisfactory correspondence. For visco us stirring the 



agreem ent is not so surprising, since we have calibrated the expressions against lOhtsuki et al 



(|2002l ). For dynamical friction, however, we have not done so, and the close match vindicates our 
geometrical model. Note, finally, that the accuracy of our geometrical model breaks down near 
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Fig. 21. — . The dimensionless stirring rates for eccentricities ((P)vs, (-P)df) and inclinations {{Q)vs, 
{Q)di) as fun ction of Hill eccentric ity, eh = v/vh- The thick gray lines represents the stirring rates 
according to lOhtsuki et al.l (j2002l ) , while the thin black lines are obtained using our geometrical 
model. The dotted line at Sh = 2.5 identifies the transition between the s.d. and d.d. -regimes. 

V ~ 2.5vh {cfi ~ 2.5) which indicates the transition between the s.d.- and d.d. -regimes, which is in 
our case sharp by construction. 



C. Calculation of interaction rates in multi-zone setting 

The interaction rates between two particle swarms over a region of space V is given by 

Ai2 = / d^x ni(x)n2(x)fTint(x)i;a, [s^^] 
Jv 



(CI) 



where ni(x) and n2(x) are the particle number densities, (Tint = t^RxRz the interaction cross section, 
Va the approach (relative) velocity, and V the volume that is under consideration. We approximate 
the ni(x) as 

Pi{x) 



ni(x) = ni{x,y,z) = Ni 



(27ra)2/ii ' 



i\^\<hi) 



(C2) 
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where A^i is the total number of particle in the swarm, hi the scaleheight, Pi{x) the probability 
density of the first particle over the radial coordinate x. The number density is zero for heights 
l^l > hi. Using Eq. IC2I we rewrite Eq. ICll as 

where the integration of the azimuthal coordinate (y-direction) has canceled one factor of Ina and 
(f)x, (jjz^ the filling factors, are defined as 

<Px = Jdx 2Rx{x)Pi{x)P2{x)- (C4) 

dz (C5) 
2hi2h2 

The filling factors have the following geometrical interpretation: it gives the fraction of the space 
with respect to radial or vertical dimension of V over which the interaction can take place. We 
therefore have that < 4>x,z < 1- If the filling factor is unity, the interaction between the two 
particle swarms can take place at any point. For example, if Rz ^ hi, /12 the particles can interact 
irrespective of their z-coordinates and </>2 = 1. However, in the more general case of fractional 
filling factors the interactions are restricted to take place in a fraction of the volume that is under 
consideration - i.e., when the particles are within a distance Rz- Let us first consider Eq. IC5I 
The integration over z proceeds over a length given by the minimum of the scaleheights, since 
interactions can only take place if both particles are present. Therefore, Eq. IC5I integrates to 
Rz/hes, where the effective scaleheight, as defined before, is the maximum of the two scaleheights. 
However, there is one caveat: (j)z is not allowed to exceed unity (indeed, we have restricted the 
interaction range before). Thus we have 

(j)z = mm{he{{/Rz, 1). (C6) 



Next, we consider the radial filling factor (p^- The situation is somewhat more complicated 
here, since the particle distributions are not centered at the same position X, see Fig. [22j In 
Fig. [22k (top panel) the distribution density function of two groups is given. The unit of length 
is arbitrary. It is assumed that the distribution is uniform with horizontal widths hxi and hx2, 
respectively. In the figure the first particle group is, without loss of generality, centered at X = 
and characterized by a scalelength hxi = 0.1. The second particle swarm is extended over a total 
length of unity and centered at X = 0.3. 

We intuitively recognize from Fig. [22] that the amount of overlap between the populations is 
given by the distance distribution P{d) among the bodies. That is, the filling factor (px represents 
the fraction of particles that comes within reach of the interaction radius -Rint- This interpreta- 
tion refines and supersedes Eq. IC41 The procedure to compute (px is outlined in Fig. [22b . the 
distribution of relative distances can be obtained from 

P{b) = [ Pi{x)P2{x + b)dx. (C7) 
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-0.2 0.0 0.2 
Position X [upper] 



0.4 0.6 0.8 
Distance b [lower] 



Fig. 22. — . (top panel) Position density distribution of two particle groups. The distributions 

are assumed to be uniform with scalelengths h^i and hx2, respectively, {bottom panel) The distance 
probability density function obtained from the cross-correlation of Pi{X) and P2{X) (black curve). 
Distributions integrate to unity. The filling factors for two interaction radii i?int,a and -Rint,b are 
indicated. 



Equation IC7I is mathematically equivalent to the cross correlation of the distributions. For the 
parameters given above, P{b) is given in the lower panel of Fig. [22] by the black solid line. The 
range in distances, or impact parameters, spans from b = —0.3 to d = 0.9, corresponding to the 
cases where the particles are furthest apart (negative values mean here that the second particle is to 
the left of the first one). P{b) also integrates to unity. The arrow in Fig. [22] indicates the interaction 
radius, Rmt- Two interaction radii are drawn. The first, Rint,a operates at distances \b\ < 0.1, while 
Pint.b does not operate at distances |6| < 0.3. If particles come within the interaction range, the 
particles are allowed to interact. In this particular example we have chosen -Rint.a = 0.05 and (px,a 
integrates to 0.1: the area of the shaded box in Fig. [22] 

In the special case where the distributions are centered at X = we have, like the ^-direction, 
P{x) = l/2hx and arrive at (px = Rx/hx,eS (where hx^eS is the largest scalelength) . Inserting 
these expressions into Eg. IC3I we obtain A12 = 7rUaiViiV2-Rx-Rz/87ra/ieff^x,eff- If we consider the 
interaction of a single protoplanet (A^i = 1) and a swarm planetesimals, the surface density of the 
latter is Ns2 = N2/ {2'Ka){2hx,cS.) and Eq. IBll is retrieved. 

The advantages of the filling factor approach is that the formalism can be extended to more 
general situations where it concerns the particle distributions or interaction radius i?int- The filling 
factor formalism becomes especially advantageous for interactions in the Hill regime, where only 
interactions at 1.7 < b/Rh < 2.5 are allowed (particles on smaller impact parameters move on 
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horseshoe orbits and do not enter the Hill sphere). The algorithm then neatly takes account of 
this inner 'gap', as is sketched in Fig. [22l In fact the function P{b) can be expressed as a series of 
step functions, that can be integrated analytically, which gives a cumulative density function, from 
which the filling factors are readily obtained. 



D. List of frequently used symbols and abbreviations 
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Table 8. List of symbols and abbreviations 



Symbol 



Description 



Reference 



a 

n 

7 

Aa 

AUint 

e 

K 

Af) 

Ps 

Pg 

Cint 
S 

^vs,50 



RB 

a 

afr 

d.d. 

e 

/vs 

/vs— 

G 
GF 

5* 



ratio R/Rh 

ratio inclinationieccentricity (= «/e = Vz/v) 

turbulent stirring parameter 

grid resolution 

total simulation width 

velocity change upon interaction 

nebula pressure parameter 

deflection angle 

accretion rate index as in dM/dt oc A4'^ 

single particle interaction rate 

group interaction rate 

internal density of bodies (including pores) 

gas density 

interaction cross section (= ttR^Rz) 
surface density in solids 
surface density in gas 

filling factor for the bodies that comprise 50% of the stirring 
power 

filling factor for interactions 
local orbital frequency 
calibration constant 
representative body 
semi-major axis 
fragment size 
impact parameter 
drag coefficient 
sound speed 
dispersion-dominated 
eccentricity (= v/vk) 
Hill eccentricity (= v/vh) 

order of unity factor that enters into the viscous stirring rate 
expression for the planar component (eccentricities) 
order of unity factor that enters into the viscous stirring rate 
expression for the vertical component (inclinations) 
Newton's constant 
gravitational focusing 

number of physical bodies belonging to species i 
effective scaleheight of interaction (= Wz/^) 
scalelength over which particles are distributed 
scaleheight 
inclination 



Eq. a 



Sect. [2X3] 
Fig.E] 
Fig. [2] 

Appendix IB. 21 
Table [2] 
Sect. [Ml 
Sect. M\ 
Eq. ED 
Eq. [CI] 
Sect. [2311 
Sect. [2311 



Sect. [2311 
Sect. [2311 
Sect. 



Appendix [C] 



Sect. [BH 
Sect. 



Sect. 



Sect. [2311 
Sect. [2311 

Sect. [211 
Appendix [BT3] 
Appendix [BXI 

Appendix [BXI 



Appendix [ 

Sect. [231 
Fig. [2] 
Fig. [2] 
Sect.[2l] 
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Table 8 — Continued 



oyniuoi 


uebci iptioii 




m 


particle mass 






dimensionless runaway growth timescale (= Tr(,/triin) 


Sect. ED 


M 


particle mass (of big bodies) 




Ml 


mass of most massive body in population 


Sect. 


M2 


mass of second-most massive body in population 


Sect. 


Mtot 


total (solid) mass in system 


Sect. 


m* 


characteristic mass of distribution 


Eq. dZl 




surface density of j-particles 




Ns{m) 


surface density spectrum 


Sect. 


9 


group size (property of RB) 


Sect. 


N* 

9 


desired group size (as given by algorithm) 


Sect. [221/ Appendix El 


Nrh 


total number of RBs (computational particles) in simulation 


Sect. 


Nres 


total number of RBs per zone (= Nrh/N^o) 


Sect. O 


N^o 


number of zones 


Fig. [2] 


fin 


fraction of particles within \b\ < 2.5Rh that enters the Hill 


Appendix |B.1.3| 




sphere 






collision probabilitv of oarticles in Hill sphere 


Appendix IRL3J 


f fra,g 


fraction of mass in fragments 


Sect. [5] 


/tot 


cumulative fraction of mass in fragments 


Sect. [5] 


P 


index of power-law mass distribution 


Sect. O 




radius corresponding to peak of Mi / m.^ 


Sect. E] 




initial planetesimal radius 




Ri(t) 


radius of most massive body at time t 


Sect. 




interaction radius for dynamical friction 


Sect. El 




interaction radius 


Sect. El 


Rf 


final radius of most massive body 




Rh 


Hill radius 


Eq. E 


-^pp 


Radius of protoplanet 




T} 

rg/oli 


transition radius between runaway growth and oligarchy in 


Eq. [21 




2-component approximation 




Rjg 


combined radius (= Ri + R2) 




RtT 


radius between runaway growth and oligarchy phase 


Eq. [25] 




interaction radius for viscous stirring 


Sect. [231 


Rvs—d 


interaction radius for viscous stirring (distant interactions) 


Sect. [221 


Rx 


geometrically constrained interaction radius in planar direc- 






tion 




Rz 


geometrically constrained interaction radius in vertical di- 






rection 




RG 


runaway growth 




s.d. 


shear-dominated 




T 


accretion (growth) timescale 


Eq. [20] 


T^2D-sd 
ac 


accretion timescale in 2D, shear-dominated regime 


Eq. [23 



- 79 - 



Table 8 — Continued 



Symbol 


Description 


Reference 


T 


runaway growth timescale 


Eq. UHl 


t 


time 




^drag 


friction (stopping) time of bodies 


Sect. [23] 


^run 


fiducial timescale that depends on initial conditions only 


Eq. [23] 


V 


planar rms-velocity dispersion 


Sect. [21] 


VM 


planar rms-velocity dispersion of heaviest particle in inter- 


Sect. [21] 




action 




Va 


approach velocity 


Sect. [21] 


Vesc 


escape velocity 


Sect. [21] 


Vh 


Hill velocity 


Eq. m 


Vk 


local orbital Keplerian velocity 




Vm 


planar rms-velocity dispersion of lightest particle in interac- 


Sect. [21] 




tion 




Vx 


maximum random velocity of the particle distribution 


Sect. [H] 


Vz 


vertical rms-velocity dispersion (usually of the large 


Sect. [21] 




body(ies)) 




W 


relative rms-velocity of two particles (excluding Keplerian 


Sect. [21] 




shear) 




Wz 


relative rms-velocity of two particles (vertical components) 


Sect. [21] 



